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Abstract 

During  the  boost  phase  of  ascent,  winds  have  a  significant  impact  on  a  launch  ve¬ 
hicle’s  angle  of  attack,  and  can  induce  large  structural  loads  on  the  vehicle.  Tradi¬ 
tional  methods  for  mitigating  these  loads  involve  measuring  the  winds  prior  to  launch 
and  designing  trajectories  to  minimize  the  vehicle  angle  of  attack  (a).  The  current 
balloon-based  method  of  collecting  wind  field  information  produces  wind  profiles  with 
significant  uncertainty  due  to  the  inherent  time  delays  associated  with  balloon  mea¬ 
surement  procedures.  Managing  the  mission  risk  caused  by  these  uncertain  wind 
measurements  has  always  been  important  to  control  system  designers.  This  thesis 
will  describe  a  novel  approach  to  managing  structural  loads  through  the  combina¬ 
tion  of  a  Light  Detection  and  Ranging  (LIDAR)  wind  sensor,  and  Model  Predictive 
Control  (MPC).  LIDAR  wind  sensors  can  provide  near  real-time  wind  measurements, 
significantly  reducing  wind  uncertainty  at  launch.  MPC  takes  full  advantage  of  this 
current  wind  information  through  a  unique  combination  of  proactive  control,  con¬ 
straint  integration  and  tuning  flexibility.  This  thesis  describes  the  development  of 
two  types  of  MPC  controllers,  as  well  as  a  baseline  controller  representative  of  cur¬ 
rent  control  methods  used  by  industry.  A  complete  description  of  Model  Predictive 
Control  theory  and  derivation  of  the  necessary  control  matrices  is  included.  The  per¬ 
formance  of  each  MPC  controller  is  compared  to  that  of  the  baseline  controller  for 
a  wide  range  of  wind  profiles  from  both  the  Eastern  and  Western  U.S.  Test  Ranges. 
Both  MPC  controllers  are  shown  to  provide  reductions  of  greater  than  50%  in  a,  Qa 
and  structural  bending  moments.  In  addition,  the  effects  of  wind  measurement  delays 
and  uncertainty  on  the  performance  of  each  controller  are  investigated. 

Thesis  Supervisor;  Frederick  W.  Boelitz 

Title:  Staff,  The  Charles  Stark  Draper  Laboratory,  Inc. 

Thesis  Advisor:  John  J.  Deyst 

Title:  Professor  of  Aeronautics  and  Astronautics,  MIT 
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Chapter  1 


Introduction 


During  the  boost  phase  of  ascent,  winds  have  a  significant  impact  on  a  laimch  vehicle’s 
angle  of  attack,  and  can  induce  large  structural  loads  on  the  vehicle.  These  structural 
loads  have  a  direct  impact  on  launch  safety,  and  are  monitored  closely  during  preflight 
launch  preparations.  Over  the  past  40  years,  launch  procedures  have  evolved  which 
attempt  to  quantify  and  minimize  the  risks  that  winds  induce  upon  launch  attempts. 
These  procedures  usually  involve  several  pre-launch  balloon  releases.  The  balloons 
are  monitored  by  radar,  which  record  their  rise  through  the  atmosphere  and  the 
associated  winds  that  are  present  at  various  altitude  levels.  Because  of  the  relative 
expense  of  these  balloons,  as  well  as  the  personnel  required  to  launch  and  track  them, 
few  are  used  for  each  launch,  resulting  in  dated  wind  field  information  which  contains 
significant  uncertainty.  The  time  between  wind  state  observation  and  vehicle  launch 
is  known  as  the  lapse  time. 

To  deal  with  this  uncertainty,  two  possible  courses  of  action  have  been  taken  in 
the  past.  The  first  is  to  directly  address  the  issue  of  wind  knowledge  certainty.  This 
is  done  by  enforcing  conservative  wind  requirements  on  the  launch  decision.  This 
approach  seeks  to  limit  mission  risk  by  ensuring  that  the  possibility  of  an  unsafe 
wind  situation  developing  within  the  lapse  time  is  below  a  certain  threshold.  Thus, 
extremely  limiting  wind  requirements  are  implemented,  which  often  cause  launch 
delays  or  cancellations  due  to  marginal  or  quickly  changing  wind  conditions. 

The  second  possible  course  of  action  is  to  implement  load  management  procedures 
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or  systems  onboard  the  vehicle,  so  that  it  can  manage  wind  changes  during  flight  in 
real-time.  Several  advances  have  been  made  in  this  area  throughout  the  past  decade. 
However,  these  changes  have  only  partially  addressed  the  issue  of  wind  uncertainty. 

This  thesis  will  address  both  of  these  methods  of  dealing  with  uncertainty  by 
proposing  and  evaluating  a  load  relief  system  which  consists  of  a  predictive  controller 
using  Model  Predictive  Control  (MPC),  working  in  conjunction  with  a  real  time  wind 
sensor,  known  as  a  Doppler  Light  Detection  and  Ranging  (LIDAR)  sensor,  which  can 
provide  accurate  and  timely  wind  knowledge  to  the  controller.  This  sensor,  which  uses 
reflected  laser  pulses  to  measure  air  movement  in  the  atmosphere  could  be  located 
on  the  vehicle  itself,  on  the  ground  near  the  launch  site,  or  on  an  aircraft  flying  over 
the  launch  site.  Although  the  technical  specifics  of  this  type  of  device  are  beyond 
the  scope  of  this  thesis,  wind  sensors  with  these  capabilities  are  currently  under 
development  in  the  private  sector,  and  will  be  available  in  a  matter  of  years. 

Using  real  time  wind  knowledge  supplied  by  the  LIDAR  wind  sensor,  the  MPC 
controller  is  able  to  anticipate  wind  disturbance  changes,  and  position  the  vehicle 
attitude  in  flight  so  as  to  minimize  the  structural  impact  of  wind  changes  on  the 
vehicle.  This  capability  significantly  expands  the  acceptable  launch  envelope  of  the 
vehicle,  resulting  in  fewer  launch  delays  and  scrubs.  In  addition,  this  technology  could 
be  incorporated  into  the  design  of  future  launch  vehicles,  possibly  allowing  structural 
mass  savings  and  equivalent  increases  in  payload  capacity  in  response  to  relaxed  load 
bearing  requirements. 

1.1  Problem  Motivation 

As  stated  previously,  wind  knowledge  uncertainty  causes  excessive  conservatism  in 
the  launch  planning  process.  This  conservatism  is  responsible  for  many  delays  and 
cancellations  during  questionable  weather  conditions.  Data  covering  launch  opera¬ 
tions  from  both  the  Eastern  Test  Range  (ETR)  at  Kennedy  Space  Center,  FL,  and  the 
Western  Test  Range  (WTR)  at  Vandenburg  AFB,  CA,  during  the  last  decade  shows 
that  nearly  half  of  all  launch  cancellations  were  due  to  high  altitude  winds.  Previous 
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industry  research  efforts  have  shown  that  launch  delays  can  cost  the  customer  an 
average  of  $500,000  per  day,  simply  for  pad  access  costs.  This  does  not  take  into 
account  the  costs  associated  with  the  delays  that  other  payloads  might  experience 
due  to  prior  launch  delays. 

Although  the  problem  of  wind  uncertainty  heis  always  existed,  the  technology 
to  implement  the  proposed  load  relief  solution  has  only  recentty  become  available. 
LIDAR  wind  sensors  with  enough  power  to  monitor  winds  up  to  altitudes  of  30 
km  are  currently  being  developed.  In  addition,  the  computing  power  required  to 
enable  real-time,  high-bandwidth  applications  of  the  computationally  intensive  MFC 
approach  have  onlj^  become  available  within  the  past  few  years.  Prior  to  that,  MFC 
was  widely  used  by  the  process  control  industry,  where  requirements  could  be  satisfied 
with  smaller  bandwidths  and  slower  cycle  times.  Continuing  growth  in  processor 
speeds  should  allow  MFC  to  be  applied  to  an  increasingly  large  segment  of  flight 
control  applications  in  the  future.  Because  of  MFC’s  predictive  nature,  it  is  able  to 
billy  utilize  current  wind  profile  information,  which  can  only  be  provided  by  a  real¬ 
time  wind  sensor,  such  as  a  Doppler  LIDAR.  Thus,  these  two  technologies  naturally 
compliment  each  other  as  an  integrated  solution  to  the  launch  uncertainty  problem. 

1.2  Overview 

The  overall  goal  of  this  research  is  to  showcase  the  advantages  and  disadvantages 
of  the  proposed  MPC/LIDAR  load-relief  system,  as  compared  to  a  traditional  load- 
relief  system  using  traditional  wind  measurement  techniques.  The  included  results 
will  demonstrate  some  of  the  significant  potential  that  MFC  and  Doppler  LIDAR  may 
have  for  contributing  to  increased  launch  probability.  In  addition,  the  MPC/LIDAR 
load-relief  system  will  be  compared  with  a  traditional  load-relief  system  coupled  with 
a  Doppler  LIDAR  wind  sensor.  This  comparison  will  highlight  the  benefits  that  MFC 
alone  can  offer  towards  increasing  launch  capability. 

The  models  used  in  these  comparisons  are  identical,  except  for  the  individual 
controller  integrated  into  each  model.  Each  vehicle  model  is  a  highly  accurate,  6 
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Degree-of-Freedom  (6DOF)  non-linear  representation  of  the  Kistler  Aerospace  Cor¬ 
poration  (KAC)  K-1  reusable  launch  vehicle.  Each  controller  guides  a  vehicle  model 
along  trajectories  generated  with  a  preflight  trajectory  generator.  Each  simulation  is 
subject  to  a  set  of  actual  wind  profiles  which  were  recorded  at  either  the  ETR  or  the 
WTR.  Over  3000  wind  profiles  were  processed  for  this  research. 

1.3  Results  Portability 

While  the  model  employed  for  this  research  represented  the  K-1  reusable  launch 
vehicle,  the  results  of  this  thesis  are  certainly  applicable  to  a  wide  range  of  launch 
vehicles,  both  reusable  and  expendable.  The  dynamics  of  most  launch  vehicles  depend 
on  a  similar  set  of  vehicle  parameters,  which  were  included  in  the  K-1  vehicle  model. 
Although  specifics  such  as  the  location  of  the  LIDAR  wind  sensor  could  vary  with 
vehicle,  the  general  result  trends  of  this  research  should  be  applicable  to  most  current 
launch  vehicle  platforms. 

1.4  Thesis  Preview 

Following  the  introductory  chapter,  this  thesis  is  organized  as  follows: 

Chapter  2  presents  the  K-1  vehicle  model  framework,  including  the  assmnptions 
and  simplifications  which  were  made  for  this  research.  Other  model  details,  such  as 
coordinate  frames,  mass  properties  and  plant  identification  are  discussed. 

Chapter  3  presents  the  mathematical  background  of  MFC  controller  theory,  as  well 
as  a  verification  of  the  MFC  controller  used  for  this  research.  Additional  features  of 
MFC  are  also  discussed,  such  as  constraints  and  slack  variables. 

Chapter  4  presents  the  architecture  development  process  which  was  used  to  choose 
the  two  MFC  control  approaches  used  in  this  research  from  an  original  set  of  five  MFC 
control  architectures. 

Chapter  5  presents  the  wind  models  used  in  this  research,  as  well  as  the  prepro¬ 
cessing  which  was  required  before  the  wind  data  could  be  used  by  the  simulations. 
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In  addition,  time  lapse  wind  processing  and  the  Trajectory  Generator  are  discussed. 

Chapter  6  presents  the  common  PID  (Proportional-plus-Integrator-plus-Derivative) 
controller  layout,  which  was  used  in  all  three  simulations  in  some  form.  The  baseline 
controller  design  process  is  discussed  in  detail. 

Chapter  7  presents  the  MPC  controller  parameter  development  process.  This 
includes  parameter  selection  methods  for  input  and  output  weighting  matrix  values, 
as  well  as  projection  and  control  horizons,  prediction  step  sizes,  controller  rates  and 
Jm  matrices. 

Chapter  8  presents  the  results  from  the  comparison  of  the  MPC  controllers  and 
the  basehne  controller  for  a  range  of  wind  profiles  and  time  lapse  scenarios. 

Chapter  9  summarizes  the  findings  of  this  research  and  concludes  with  lessons 
learned  and  recommendations  for  future  research. 

A  comprehensive  list  of  symbols,  acronyms  and  conventions  is  included  in  Ap¬ 
pendix  A. 
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Chapter  2 


Vehicle  Model 

2.1  Model  Framework 

Traditionally,  vehicle  simulations  are  constructed  in  either  the  FORTRAN,  C,  or 
C++  software  languages  as  collections  of  stand-alone  executable  files.  Although  this 
method  produces  simulations  with  fast  execution  times,  quickly  modifying  or  ad¬ 
justing  the  simulations  becomes  increasingly  difficult  with  increased  size.  The  cost 
required  to  develop  the  input/output  infrastructure  needed  to  support  rapid  engi¬ 
neering  analysis  and  design  is  substantial.  One  advantage  of  stand-alone  simulations, 
however,  is  the  supply  of  reliable  functions  that  have  been  previously  coded  and  tested 
in  other  programs,  such  as  standard  atmosphere  models  or  gravity  models. 

The  standard  modelling  and  simulation  tool  employed  today  is  SIMULINK,  by 
Math  Works,  Inc.  In  addition  to  providing  a  graphically  driven  user  interface  with 
strong  input /output  support,  the  software  package  can  also  work  in  unison  with 
Math  Works’  other  main  software  product,  MATLAB.  MATLAB  provides  basic  lin¬ 
ear  algebra  functionality  along  with  other  tools  for  designing  control  systems  and 
displaying  signals.  The  tools  are  provided  in  the  form  of  software  toolboxes  that  can 
be  purchased  separately  from  MathWorks. 

The  ideal  simulation  environment  would  couple  the  rapid  simulation  development 
capabilities  supplied  by  Simulink,  with  the  speed  supplied  by  stand-alone  C  based 
functions.  Fortunately,  Simulink  and  MATLAB  provide  an  open  software  interface 
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that  allows  custom  designed  code  to  be  executed  from  within  its  environment. 

A  Simulink  S-Function  (system-function)  is  a  computer  language  description  of  a 
Simulink  function  block.  S-Functions  can  be  written  in  MATLAB,  C,  C-|— f,  Ada,  or 
FORTRAN.  S-Functions  use  a  special  calling  syntax  that  enables  custom  designed 
functions  to  interact  with  Simulink’s  equation  solvers.  Typically,  a  programmer  cre¬ 
ates  a  wrapper  function  that  calls  a  separate  model  function.  The  wrapper  function 
contains  all  of  the  required  application  programming  interface  functions  required  by 
the  Simulink  environment,  while  the  model  function  contains  the  actual  executable 
code.  In  this  way,  the  model  function  (e.g.  1976  US  Standard  Atmosphere)  is  neatly 
isolated  from  Simulink,  thereby  making  it  easily  portable  to  other  simulation  envi¬ 
ronments. 


2.2  Non-Linear  Vehicle  Model 

For  this  thesis,  all  of  the  dynamic  models  that  collectively  describe  the  motion  of  the 
Kistler  K-1  are  coded  as  separate  C  S-Functions.  These  functions  are  integrated  into 
Simulink  by  calling  them  from  their  corresponding  wrapper  functions,  identified  by 
the  naming  syntax  filename. SMLK.c.  Hereafter,  all  S-Functions  will  be  referred 
to  simply  by  their  filenames.  All  of  the  S-Functions  are  compiled  into  dynamic  link 
libraries  {.dll)  on  a  PC  using  Microsoft’s  Visual  C-f- 1-  Studio  version  6.0.  The  vector 
and  transformation  matrix  naming  conventions  used  in  this  thesis  are  detailed  in 
Appendix  A.  A  more  complete  list  of  the  equations  contained  in  each  S-Function  can 
be  found  in  Appendix  B.  The  critical  dynamic  models  used  to  describe  the  motion 
of  any  rigid  body  vehicle  travelling  about  a  spherical  earth  include: 

•  accel  :  Computes  the  r;  as  well  as  the  The  primary  inputs  to  the 

model  are  the  vehicles  mass  properties  and  the  forces  and  torques  acting  on  the 
body. 

•  eulerRates  :  Computes  the  time  rate  of  change  of  the  yaw  (^/;) ,  pitch  (0) ,  and 

roll  (0)  euler  angle  sequence.  The  input  is  the  computed  by  the  accel 
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function  and  the  current  9,  0,  euler  angle  sequence. 


•  envGMST  :  Computes  transformation  matrices  Tf  ^e-  inputs  to  this 
model  are  time  and  the  rotation  rate  of  the  earth. 

•  xformEuler  :  Computes  transformation  matrices  T^,  T^,  Tf^,  ^ 

and  Ti'g.  These  transformations  are  computed  using  the  ciurrent  vehicle  latitude, 
longitude,  attitude  and  Tg. 

•  velocityOmega  ;  Computes  various  velocities  and  angular  rates  of  the  body 

with  respect  to  a  variety  of  frames.  These  include  ^e|g,  Ib’ 

inputs  to  the  model  include  the 
vehicle  longitude,  T§,  tI^,  '^l\p  ~pI\j,  rotation  rate 

of  the  earth. 

•  gravityShepperd  :  Computes  the  force  vector  on  the  body  due  to  gravity, 
expressed  in  the  body  frame,  F the  gravity  vector  in  both  the  earth  and 
body  frames,  ~Vg\^  and  ~Vg\^  respectively,  and  the  weight  of  the  vehicle.  The 
inputs  to  the  model  are  the  vehicle  mass,  the  and  the 

•  aeroParameters  ;  Computes  the  angle  of  attack  (a),  sideslip  angle  (/?),  total 

angle  of  attack  (o;^),  aerodynamic  roll  angle  (0*),  airspeed,  dynamic  pressure 
((5),  and  mach  number  based  on  the  local  air  density  (p),  and  speed  of 

sound  (SoS). 

•  atmosphereUS76  :  Computes  the  local  pressure,  temperature,  density,  and 
speed  of  sound  as  a  function  of  the  local  radius  of  the  Earth  and  the  altitude  of 
the  vehicle  above  mean  sea  level. 

•  xyz21at_lon_2ilt  :  Computes  the  latitude,  longitude,  altitude  and  the  local 

radius  of  the  earth,  re,  as  a  function  of  the  Parameters  include  the 

equatorial  earth  radius,  reO,  and  the  earth  flattening  constant,  k. 

In  order  to  uniquely  deflne  the  K-1,  the  following  additional  models  are  required: 
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•  engine  :  Computes  the  force  and  torque  vectors  produced  by  the  rocket  engines 
expressed  in  the  body  frame.  The  inputs  to  the  model  are  the  location  of  the 
engine  hinge  point  in  the  body  frame,  the  local  air  pressure,  the  current  throttle 
setting  [0-1],  the  vehicle  CG  location,  and  the  pitch  and  yaw  gimbal  angles  with 
respect  to  the  bod)"  frame,  5p  and  5y. 

•  aeroKlBoost  :  Computes  the  aerodj^namic  force  and  torque  vectors  based 
upon  the  current  Q,  (p*,  CG  location,  and  mach  number. 

The  full  non-linear  K-1  simulation  is  comprised  of  approximately  fifteen  S-Functions 
linked  together  to  form  a  comprehensive  6  Degree-of-Freedom  rigid  body  simulation. 
A  collection  of  pictures  and  descriptions  of  all  6DOF  models  used  in  this  research 
is  included  in  Appendix  E.  Due  to  the  flexibility  of  the  Simulink  environment,  any 
number  of  continuous  integrators  can  be  chosen.  (r/c4  was  typically  used  with  a  sam¬ 
ple  time  of  At  =  0.02  sec)  depending  upon  the  speed  and  fidelity  required.  All  work 
described  in  this  thesis  was  completed  using  Simulink  version  4  and  MATLAB  version 
6.1.0.450  release  12.1. 

2.2.1  Standard  Coordinate  Frame  Descriptions 

The  full  simulation  uses  the  four  standard  frames  found  in  space  system  simulations. 
All  frames  are  orthonormal  and  right-handed.  Figures  2-1  and  2-2  illustrate  these 
coordinate  frames. 

•  The  earth-centered  inertial  frame  (ECI),  commonly  called  the  Inertial  Frame. 
This  frame,  with  axes  7,  J  and  /v,  is  fixed  at  the  earth’s  center,  with  the  I 
axis  pointing  towards  the  vernal  equinox,  the  J  axis  perpendicular  to  the  I  axis 
within  the  equatorial  plane,  and  the  K  axis  completing  the  right-hand  system, 
pointing  out  the  North  Pole,  along  the  earth’s  axis  of  rotation. 

•  The  earth-centered  rotating  frame  (ECR),  commonly  called  the  Earth  Frame. 
This  frame  is  similar  to  the  ECI  frame,  except  that  it  rotates  with  the  earth. 
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with  the  primary  axis  aligned  with  a  fixed  meridian,  normally  the  Greenwich 
Meridian. 

•  The  local  geographic  frame,  commonly  called  the  NED  or  LG  Frame.  This 
frame  is  fixed  at  some  point  on  the  earth’s  surface,  with  the  primary  axis 
pointing  North,  the  second  pointing  East  and  the  third  axis,  completing  the 
right-handed  system,  pointing  Down. 

•  The  body-fixed  frame,  commonly  called  the  Body  Frame.  This  frame  is  fixed 
at  some  point  on  the  body,  usually  the  CG.  The  primary  axis  (X)  points  out 
the  nose  of  the  vehicle,  with  the  second  axis  {Y)  perpendicular  to  the  primary 
axis,  usually  out  the  right  side  of  the  vehicle.  The  third  axis  (Z)  completes  the 
right-handed  system,  normally  pointing  out  the  bottom  of  the  vehicle. 


1^ 

*^ECI,ECR 


Figure  2-1;  ECI,  ECR  and  NED  Coordinate  Frames 
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Figure  2-2:  Body  Coordinate  Frame 


2.2.2  Assumptions  and  Simplifications 

Several  assumptions  and  simplifications  were  made  in  constructing  the  non-linear 
vehicle  model.  These  are  listed  below  and  on  the  following  pages,  along  with  their 
justifications: 

•  The  Earth  is  modelled  as  a  flat  surface  with  an  X,  Y,  Z  inertial  reference  frame. 
Only  this  frame  and  the  body  frame  are  required.  Because  this  research  is 
constrained  to  the  relatively  short  (134  sec)  boost  phase  of  ascent,  the  vehicle 
travels  over  a  very  small  percentage  of  the  earth’s  surface.  Therefore,  any  errors 
resulting  from  differences  in  navigating  through  a  spherical  Earth  environment 
as  opposed  to  a  flat  Earth  environment  can  be  shown  to  be  negligible. 

•  The  vehicle  is  assumed  to  be  a  rigid  body.  Effects  due  to  liquid  fuel  slosh 
and  structural  bending  are  assumed  to  be  higher  order  and  will  not  impact 
the  performance  of  the  proposed  guidance  system.  This  assumption  is  based 
on  previous  analysis  of  the  K-1  which  shows  that  the  first  structural  bending 
mode  occurs  at  approximately  18.5  rad/sec,  and  the  first  slosh  mode  occurs 
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at  approximately  2.3  rad/sec  [3].  To  capture  the  real  world  effects  that  these 
properties  would  have  on  system  bandwidth,  the  control  system  design  for  this 
research  was  limited  to  between  1  and  2  rad/sec,  a  typical  range  for  a  launch 
vehicle. 

•  Engine  gimbal  dynamics  are  assumed  to  be  hnear  and  second  order  in  nature. 
Frequency  domain  design  specifications  for  the  K-1  rocket  gimbal  system  show 
that  this  assumption  is  sufficient  to  model  the  system’s  expected  response. 

•  Tail-Wags-Dog  (TWD)  effects  are  excluded  from  the  simulation.  TWD  refers 
to  the  dynamic  effects  due  to  the  movement  of  the  main  engine  nozzle  which  is 
attached  to  the  aft  end  of  the  vehicle  body.  Previous  analysis  of  the  K-1  has 
shown  that  the  frequency  of  the  TWD  zero  is  much  higher  than  the  bandwidth 
of  the  proposed  control  system,  at  approximately  30  rad/sec  [3].  As  with  body 
bending,  this  is  assumed  to  be  a  higher  order  effect  that  will  not  influence  the 
performance  of  the  proposed  load  rehef  system. 

•  Any  transportation  lags  associated  with  on-board  processing  are  assumed  to  be 
negligible.  For  the  K-1,  this  is  not  necessarily  true.  However,  the  effect  that 
computation  lags  may  have  on  system  performance  is  acknowledged  by  limiting 
the  bandwidth  of  the  proposed  control  system  as  previously  stated. 

•  The  K-1  has  three  Russian  NK-33  booster  engines,  each  delivering  approxi¬ 
mately  395,000  Ibf  of  thrust,  when  measured  in  the  vacuum  of  space.  Pitch 
and  yaw  control  is  achieved  by  simultaneously  gimballing  all  three  engines  in 
the  same  direction.  Roll  control  is  achieved  by  differentially  gimballing  the 
outboard  engines.  For  this  research,  all  three  engines  are  represented  as  a  sin¬ 
gle,  central  engine  of  equivalent  thrust.  Pitch  and  yaw  control  methods  remain 
unchanged.  However,  because  the  simulation  vehicle  has  only  one  engine,  roll 
control  is  not  possible.  In  all  other  respects,  a  single  engine  can  effectively 
simulate  the  same  dynamic  effect  as  three  separate  engines. 

•  Although  the  real  K-1  has  aerodynamic  damping  forces,  these  are  not  modelled 
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in  the  simulation.  For  the  K-1,  these  forces  are  small  and  their  tendency  is 
to  add  robustness  to  the  control  design.  Eliminating  this  effect  adds  a  small 
measure  of  conservatism  to  the  analysis. 

•  The  actual  K-1  cross  products  of  inertia,  while  not  exactly  zero,  are  so  small 
that  they  can  safely  be  assumed  to  be  equal  to  zero  without  impacting  the 
validity  of  the  vehicle  simulation.  The  average  I^y  inertia,  for  example,  is  only 
2.7%  of  the  average  roll  inertia,  Ixx  0.03%  of  the  average  pitch  inertia,  lyy. 

•  The  center  of  gravity  of  the  body  is  assumed  to  move  as  a  function  of  the  vehicle 
mass,  and  its  motion  is  constrained  to  the  body  x-axis.  This  assumption  is  valid 
because  the  actual  off-axis  component  of  the  CG  location  is  extremely  small, 
on  the  order  of  1  inch  for  a  vehicle  which  is  approximately  115  feet  in  length 
and  22  feet  in  diameter. 

•  Feedback  of  angular  rate,  attitude,  and  position  into  the  control  system  is  as¬ 
sumed  to  be  perfect.  Although  all  measurement  devices  produce  biased  and 
noisy  signals,  the  issues  caused  by  these  inaccuracies  are  beyond  the  scope  of 
this  thesis.  In  addition,  because  all  systems  will  be  affected  in  similar  ways 
by  signal  noise,  the  comparison  framework  of  this  thesis  will  act  to  minimize 
the  impact  of  this  noise  on  the  relative  performance  results  of  the  baseline  and 
proposed  load  relief  guidance  systems. 

•  All  trajectories  used  in  this  research  are  solely  pitching  trajectories,  and  include 
no  non-zero  yaw  or  roll  trajectory  components.  This  allows  traditional  control 
methods  to  be  used  for  yaw  plane  control,  and  negates  the  need  for  roll  control. 

•  All  winds  are  along  the  inertial  frame  Z  axis,  as  shown  in  Figure  2-3.  There  are 
no  lateral  or  vertical  components  to  wind  disturbances. 

•  Because  wind  disturbances  only  occur  in  the  pitch  plane,  all  research  and  anal¬ 
ysis  was  conducted  solely  in  this  plane. 

•  Winds  are  assumed  to  be  constant  at  any  position  at  a  given  altitude. 
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•  All  simulation  wind  profiles  are  assumed  to  accurately  represent  actual  wind 
profiles,  although  they  were  gathered  from  balloon  measurements,  which  have 
been  shown  to  have  a  slight  dampening  effect  on  wind  measurements. 


2.3  Simplified  Flat-Earth  Non-Linear  Vehicle  Model 

As  previously  stated,  to  reduce  complexity,  increase  simulation  speed  and  provide 
for  an  easier  to  understand  and  operate  framework,  the  simulation  environment 
was  converted  from  a  Round-Earth  to  a  Flat-Earth  environment.  Only  two  coor¬ 
dinate  frames  were  retained,  as  shown  in  Figure  2-3.  The  S-Functions  envGMST 
and  xyz2latJon.alt  were  removed.  In  addition,  the  S-Functions  xformEuler  and 
velodtyOmega  were  rewritten  to  eliminate  all  unnecessary  calculations  resulting  from 
the  simulation  environment  transition. 

2.3.1  Coordinate  Frame  Descriptions 

The  simplified  vehicle  simulation  contains  two  frames,  the  inertial  frame  and  the  body 
frame.  The  inertial  frame  is  fixed  at  the  launch  site  and  follows  standard  right-hand, 
orthonormal  conventions.  The  X  axis  points  up,  and  is  perpendicular  to  the  Y  and 
Z  axes,  which  lie  in  the  plane  of  the  Earth’s  surface,  following  the  right  hand  rule. 
The  second  frame  is  the  body  frame,  which  is  fixed  to  the  body  with  its  origin  at  the 
vehicle’s  CG,  in  the  traditional  nose,  right  wing,  down  configuration  as  described 
previously.  Thus,  the  two  frames  can  begin  the  flight  co-located  and  aligned.  A 
picture  of  the  inertial  and  body  frames  is  included  in  Figure  2-3.  Figure  2-4  describes 
the  variables,  angle  of  attack  (a),  pitch  (0)  and  flight  path  angle  (7).  represents 
the  velocity  vector  of  the  body,  with  respect  to  the  air.  This  notation  convention  is 
explained  in  more  detail  in  Appendix  A. 


33 


2.3.2  Vehicle  Mass  Properties 

The  vehicle  inertia  tensor  and  CG  location  information  used  in  this  simulation  are 
representative  of  a  Kistler  K-1.  Table  2.1  contains  the  actual  values  used  in  the  ve¬ 
hicle  simulation.  Although  not  shown,  the  magnitudes  of  the  oflF-axis  inertia  tensor 
components  and  the  Y  and  Z  CG  positions  are  quite  small,  as  stated  in  the  Assump¬ 
tions  and  Simplifications  section.  In  addition,  because  all  analysis  is  done  in  the  pitch 
plane,  and  the  vehicle  has  virtually  no  roll  or  yaw  inputs  during  flight,  off-axis  inertia 
components  would  not  significantly  affect  performance,  regardless  of  their  size.  The 
CG  position  is  measured  from  a  point  fixed  1000  inches  behind  the  vehicle  separation 
plane,  where  the  two  stages  of  the  K-1  connect. 


Mass  (slugs) 

CGx  (ft) 

Ixx  (slug-ft^) 

lyy 

Izz 

26243 

73.53 

212,407 

19,183,122 

19,192,511 

90%  Fueled 

24825 

74.17 

18,939,892 

18,949,259 

80%  Fueled 

23407 

74.98 

18,602,901 

18,612,290 

70%  Fueled 

21989 

75.98 

201,450 

18,171,523 

18,180,934 

60%  Fueled 

20571 

17,631,988 

17,641,377 

50%  Fueled 

19153 

194,140 

16,960,056 

16,969,467 

16,123,094 

16,132,526 

■gifeiiaireiail 

186,821 

15,074,392 

15,083,824 

20%  Fueled 

14899 

85.33 

183,158 

13,746,394 

13,755,869 

10%  Fueled 

13481 

88.74 

179,490 

12,034,330 

12,043,805 

Table  2.1:  Vehicle  Simulation  Mass  Properties 


2.3.3  Plant  Creation 

A  separate  plant  ID  program  was  used  to  create  a  linear  model  of  the  vehicle,  updated 
at  preset  increments  throughout  a  nominal,  no-wind  flight.  This  linear  representation 
of  the  vehicle  pitch  dynamics  remained  valid  in  the  face  of  small  variations  in  pitch 
or  angle-of-attack  that  might  exist  throughout  any  given  ascent  trajectory  due  to 
wind  variations.  This  data  was  saved  in  the  form  of  the  finear,  time-invariant  (LTI), 
discrete  A,  B,  C  and  D  state-spane  matrices  shown  in  Equation  3.1. 
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Chapter  3 


MPC  Controller  Theory 

3.1  MPC  Overview 

Model  Predictive  Control,  as  a  general  concept,  has  existed  for  many  years.  Because  of 
MFC’s  large  computational  requirements,  widespread  interest  first  began  during  the 
1980’s  in  applications  with  low-rate  control  requirements  and  low-bandwidths,  such  as 
those  processes  found  in  the  chemical  and  process  control  industries  [2] .  Only  recently, 
with  significant  increases  in  computer  processor  speeds,  has  it  become  feasible  to 
apply  MPC  to  high  bandwidth  systems  such  as  flight  vehicles,  which  require  high- 
rate  control  processes. 


d(t) 


Figure  3-1:  MPC  Operation  Graphical  Overview 
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The  conceptual  structure  of  MFC  is  shown  in  Figure  3-1.  As  its  name  implies, 
MFC  uses  an  internal  model  and  knowledge  of  future  and  present  measured  and  esti¬ 
mated  unmeasured  disturbances,  reference  trajectories  and  current  states  to  predict 
the  output  of  a  given  system  over  a  finite  time  horizon,  known  as  the  prediction, 
or  output  horizon,  p.  Recent  wind  studies  have  shown  that  wind  behavior  does  not 
change  in  any  statistically  significant  way  over  relatively  short  times,  such  as  five 
minutes  or  less.  Since  current  LIDAR  wind  sensors  are  capable  of  measuring  the 
wind  field  in  less  than  five  minutes,  they  can  provide  “future”  disturbance  informa¬ 
tion  to  a  MFC  controller  by  simply  providing  the  current  wind  information  at  the 
location  where  the  vehicle  will  be  at  the  appropriate  times  in  the  future.  Because  this 
research  utilizes  full  disturbance  information  provided  by  the  LIDAR  wind  sensor. 
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all  unmeasured  disturbance  estimates  are  neglected.  However,  for  completeness  sake, 
the  unmeasured  disturbance  terms  are  included  in  the  MFC  matrix  derivations.  A 
graphical  description  of  this  process  is  shown  in  Figure  3-2.  This  prediction  ability 
allows  MFC  to  solve  an  optimal  control  problem  on-line  in  order  to  minimize  the 
error  between  the  predicted  output  and  a  desired  output  trajectory,  also  known  as 
the  reference  trajectory.  This  process  could  be  subject  to  constraints  on  the  ma¬ 
nipulated  inputs  and  outputs.  The  result  of  the  optimization  is  a  series  of  optimal 
input  commands,  each  applying  to  a  given  time  step  within  a  second  horizon,  called 
the  control  or  input  horizon,  m,  also  shown  in  Figure  3-2.  Between  the  control  and 
prediction  horizon,  the  input  command  is  held  constant.  At  each  time  step,  only  the 
first  input  command  is  applied  to  the  system.  The  remainder  of  the  optimal  input 
commands  are  discarded  and  the  system  is  propagated  forward  one  time  step.  At 
this  point  the  optimal  control  problem  is  solved  again,  with  an  appropriate  update  to 
all  measurements  from  the  plant,  as  well  as  the  applicable  disturbance  and  reference 
trajectories.  As  this  process  is  repeated  and  the  system  progresses  through  time,  the 
projection  and  control  horizons  recede  at  the  same  rate. 

3.2  Unconstrained  MPC  Optimal  Solution 

Consider  the  discrete  LTI  plant  with  the  state  and  output  equations  as  follows: 


x{k  -h  1)  =  Ax{k)  -h  Buu{k)  -h  B^v{k)  +  Bdd{k)  (3.1) 

y{k)  =  Cx{k)  +  Dvv{k)  +  Ddd{k) 

where  x(k)  e  u{k)  G  y{k)  G  v{k)  G  and  d{k)  G 

As  stated  previously,  v{k)  and  d{k)  represent  the  measured  and  estimated  unmea¬ 
sured  disturbances  experienced  by  the  plant,  respectively.  In  the  same  way,  u{k)  and 
y{k)  represent  the  input  to  and  output  from  the  plant  model  and  xik)  represents  the 
present  model  states.  Model  Fredictive  Control  minimizes  the  cost  function  shown  in 
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Equation  3.5  to  find  the  optimal  input  command  series,  in  the  presence  of  the  three 
error  terms  shown  in  that  equation  and  their  respective  weighting  matrices,  Wu,  wa 
and  Wy.  These  equations  are  shown  below. 


p-i 

^  ][u{k  +  i)  —  (3-2) 

i=0 

m— 1 

'Y^z{k  +  ifwA{i)  (3-3) 

Y^[y{k  +  i)-  r{k  +  i)fwy{i)  (3.4) 

i=0 

This  cost  function  penalizes  the  error  between  the  commanded  output  and  the 
predicted  output  (Equation  3.4),  the  trim  input  and  the  predicted  input  (Equation 
3.2),  and  the  difference  between  the  current  input  command  and  the  previous  input 
command  (Equation  3.3)  as  shown  by  Equation  3.6.  Thus,  it  exhibits  features  of  cost 
functions  found  in  both  [2]  and  [4],  regulating  both  inputs  and  differences  between 
consecutive  inputs,  z  is  called  the  optimization  parameter.  Because  the  plant  in 
this  research  only  has  one  output,  an  overline  represents  a  variable  which  is  in  vector 
form.  As  shown  in  Equation  3.5,  the  optimal  z{k)  must  be  found  to  minimize  the 
cost  function,  J.  The  optimal  ^(/c)  will  be  denoted  as  'z*{k). 


p-i 


m— 1 


J(z(k))  =  min  +  ^)  “  '^rik  +  i)]^u>„(z)  +  ^  z{k  +  2)"Wa(z)+ 


p-i 


'^[yik  +  i)  -  r{k  +  i)fwy(i)  (3.5) 


i=0 


Several  other  parameters  must  be  defined  before  continuing.  These  include  nu  and 
ny,  which  represent  the  number  of  inputs  to  and  outputs  from  the  plant,  respectively. 
nd  represents  the  number  of  estimated  unmeasured  disturbances  and  nv  represents 
the  number  of  measured  disturbances.  Finally,  nx  represents  the  number  of  model 
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states. 

The  optimization  parameter  z{k)  is  related  to  the  variation  of  the  input  variables 
as  shown  below: 


Au{k)  =  JM^k)  where  ^  ^(p*nu)x{m*nu)  (3  g) 

Au{k)  =  u{k)  —  u{k  —  1) 

Jm  is  a  matrix  used  to  impose  additional  constraints  on  the  optimum  Au{k).  The 
function  and  structure  of  this  matrix  are  discussed  in  Section  3.4,  Blocking  and  Basis 
Functions. 

To  implement  this  cost  function  in  an  effective  and  flexible  manner,  it  should  be 
rewritten  in  matrix  form.  To  do  this,  several  new  terms  must  be  defined,  as  shown 
below. 


u{k)  =  Ipu{k  —  1)  +  KiAu{k) 


(3.7) 


Ii  =  Identity{nu  x  nu) 


(3.8) 


Ip  = 


h 

h 


Ki 


/i  0  •  •  •  0 

h  h  ■■■  0 

h  h  •••  h 


(3.9) 


(3.10) 


We  can  now  write  this  cost  function  in  matrix  form: 
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J  —  min  [w(fc)  —  UT{k)]'^W^[u(k)  —  Urik)]  +  z^{k)JjfWj\JMz{k)+ 

Rl-)-f(fc)]’'H;|5(t)-r(S,-)l  (3.11) 

If  no  constraints  exist,  a  closed-form  analytical  minimum  of  the  cost  function  can 
be  found.  This  process  is  described  in  the  following  section  and  was  derived  from  an 
initial  treatment  in  [1]. 

Let  y{k+i\k)  and  Au{k+i\k)  be  the  output  and  input  change  predictions  obtained 
by  iterating  the  model  i  times  into  the  future,  from  the  current  state,  k.  As  previously 
defined,  p  and  m  are  the  prediction  and  control  horizons,  respectively.  The  vectors 
y{k)  and  Au{k)  are  structured  as  shown  below. 

y{k  +  l|/ir) 

y{k)=  +  ^^{p.ny)xl  (3J2) 

_  y{k  +  p\k)  _ 

Au{k\k) 

_  Au(/c-fl|A:)  c  ivi 

Au{k)=  (3  13) 

Au{k  +  p  —  1|A’) 

Using  this  convention,  we  can  also  write  the  prediction  of  the  input  in  vector  form, 
as  shown  below. 

u{k) 

U{k)=  +  ^3^(p*nn)xl  (3  14) 

u{k  -hp  —  1) 

The  known  input,  measured  disturbance  and  estimated  unmeasured  disturbance 
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trajectories  are  now  defined  as  shown  below. 


urik) 


UT{k) 

Urik  +  1) 

uxik  +  p-  1) 


(3.15) 


v{k)  = 


v{k) 
v(k  +  1) 

v{k+p) 


Xl 


(3.16) 


d{k)  = 


d{k) 

d{k  + 1) 


(3.17) 


[  d{k  +  p)  J 

At  this  point  we  are  ready  to  begin  calculating  y(k).  We  will  begin  with  the 
standard  discrete  LTI  state-space  equations,  as  shown  in  Equation  3.1.  This  equation 
shows  the  system  at  time  step  k. 


x(k  +  1)  =  Ax{k)  +  Buu{k)  +  Byv{k)  -h  Bdd{k)  (3.18) 

y{k)  =  Cx{k)  +  D^v{k)  +  Ddd{k) 


The  system  equations  at  time  step  k  are  written  below. 


x{k  -h  2)  =  Ax{k  +  1)  +  Byu{k  +  1)  +  B.uv{k  +  1)  +  Bdd{k  +  1)  (3.19) 

y{k  +  1)  =  Cx{k  +  1)  -1-  D^v{k  +  1)  +  Ddd{k  +  1) 


To  reach  the  complete  realization  of  these  equations,  the  terms  of  a; (A;  -|- 1)  from 
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Equation  3.18,  must  be  inserted  where  appropriate  into  Equation  3.19,  as  shown 
below. 


x{k  +  2)  =  A[Ax{k)  +  B^v{k)  +  B^v{k)  +  Bdd{k)]+ 

Bvu{k  +  1)  +  Byv{k  +  1)  +  Bdd{k  +  1) 
y{k  +  1)  =  C[Ax{k)  +  Byu{k)  +  B^vik)  +  Bdd{k)]+ 

DyV{k  +  1)  +  Dddik  +  1) 

Writing  the  system  equations  yet  again,  for  time  step  fc  +  2,  gives  the  following 


equations. 


x{k  +  3)  =  Ax{k  +  2)  +  B^u{k  +  2)  +  Bj,v{k  +  2)  +  Bdd{k  +  2)  (3.21) 

y{k  +  2)  =  Cx{k  +  2)  +  Dvv{k  +  2)  +  Ddd{k  +  2) 

As  before,  the  terms  of  x(A;  +  2)  must  be  inserted  where  appropriate  into  Equation 
3.21,  as  shown  below. 


x{k  +  3)  =  A[A[Ax{k)  +  B^,u{k)  +  Byv{k)  +  Bdd{k)]+ 

Buu{k  +  1)  +  Byv{k  +  1)  +  Bdd{k  +  1)]+  (3.22) 

Byu{k  +  2)  +  Byv(k  +  2)  +  Bdd{k  +  2) 
y(^k  +  2)  =  (^[^[^^(A;)  +  Buuiji)  +  ByV{k^  +  Bdd(^k)]A 

Byu{k  4- 1)  +  Byv{k  +  1)  +  Bdd{k  +  1)]+ 

Dyv{k  +  2)  +  Ddd{k  +  2) 

Noticing  the  emerging  pattern,  the  prediction  of  y,  at  time  k,  to  any  time  k  +  i, 
can  be  expressed  as: 


44 


y{k  +  i\k)  = 


C 


i—1 


h=0 


u{k  —  l)  “1“  ^  ^  ^u{k  “h  J )  “1“  Byv{k  /i)  “1“  B(id(^k  ■!“  /i) 
j=o 

”1“  Djjv(^k  “1“  i)  ~l“  Dnd{k  -\-  i)  (3.23) 


To  realize  this  equation  in  vector  form,  several  additional  matrices  must  be  defined, 
as  shown  below. 


y{k)  = 


y{k  +  l\k) 
y{k  +  2\k) 

y{k+p\k) 


xl 


CA 

CA^ 

CAP 


^  ^(p*n3/)> 


CB, 

CB^  +  CAB^ 


(3.24) 


(3.25) 


(3.26) 


CBu  0 

CB^  +  CABu  CBu 

EVoCA^Bn  ElllCA^B, 


0  0 
0  0 
■■■  0 
•••  CBu 


^  ^(.p*ny)x{p*nu)  (3.27) 
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0 


H^  = 


CB, 

Dy 

0 

CAB^ 

CBy 

Dy 

CAP-^B^ 

CAP-^By 

CAP-'^B, 

^  ^{p*ny)>c{p*nu+nv)  (3.28) 


Hd 


CBd 

Dd 

0 

CABd 

CBd 

Dd 

CAP-^Bd 

CAP-^Bd 

CAP-^B, 

0 

0 


G  (3.29) 


Using  these  matrices,  the  prediction  of  y{k)  can  now  be  rewritten  as  shown  below: 


y{k)  =  Sxx{k)  +  Suiu{k  —  1)  +  Su^u{k)  +  Hyv{k)  +  Hdd{k)  (3.30) 
Likewise,  the  cost  function  in  its  most  basic  matrix  form  can  now  be  written  as: 


J  =  |«(i)  -  v^(k)YW„[ti(k)  -  k?(A-)1  +  Aii’'(/t)HiAt/(A:)+ 

\y(k)  -  f(J:)]^H',|5(fc)  -  r(k)\  (3.31) 

To  find  the  analytical  minimum  of  this  cost  function,  the  derivative  of  the  function 
must  be  taken.  The  cost  function  will  be  broken  into  three  separate  parts,  Ji,  J2  and 
J3  for  simplicity.  In  addition,  the  {k)  convention  will  be  neglected  for  brevity’s  sake. 
Taking  the  first  term  of  Equation  3.31  to  be  Ji  we  have: 

Ji  =  [«  -  [u  -  ut]  (3.32) 

Substituting  from  Equation  3.7  yields: 
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Ji  =  [Ipu{—1)  +  KiAu  —  UT]^Wu[Ipu{—l)  +  KiAu  —  Ut]  (3.33) 

Substituting  once  again  from  Equation  3.6  yields: 

+  ^iJm^  ~  +  KiJmz  ~  Wx]  (3.34) 

Taking  the  partial  now  yields: 

B  7 

=  2[4u(-l)  -  utYWuK^Jm  +  2z'^JIK^W^KiJm  (3.35) 

uZ 

Taking  the  second  term  of  Equation  3.31  as  J2  and  making  the  appropriate  substitu¬ 
tion  from  Equation  3.6,  we  have: 

J2  =  =  z'^JIjWaJmz  (3.36) 

Taking  the  partial  yields: 

^  =  2z^JI,W^Jm  (3.37) 

Finally,  taking  the  third  term  of  Equation  3.31  as  J3  we  have: 

^3  =  [^  -  rfWy  [y  -  f]  (3.38) 

Inserting  Equation  3.30  where  appropriate  yields: 

J3  =  ['S'xX  -I-  5«iu(-l)  -I-  SuAu  +  Hj;V  -h  Had  -  f^Wy 

[5a;a:  +  Suiu{—1)  +  SuAu  -I-  HyV  -h  Had  -  r]  (3.39) 

Combining  all  constants,  we  can  define  a  new  term  F  as: 

F  =  [iSjja:  -f-  Su\u{~Y  4“  dlyV  -|-  Hfjd  —  r]  (3.40) 
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Inserting  F,  J3  is  now: 


=  |F  +  S.AtifWyF  +  S.Am] 


(3.41) 


Taking  the  partial  yields: 


dz 


=  2F'^WySnJM  +  2-fjljSlWySJ, 


M 


(3.42) 


Setting  the  derivatives  equal  to  zero  and  solving  will  yield  the  function  minimum. 


dJ  dJi  dJ2  dJz 
^  =  ai 


From  above,  this  is  equal  to: 


(3.43) 


[Ipu{—1)  —  ut]^WuKiJm  +  z*^  JJjK^WuKiJm  + 

F'^WySuJM  +  r'^JhSlWyS^JM  =  0  (3.44) 

Consolidating  the  1*  terms  yields: 


z*^(JJ,KIW^KiJm  +  JJjW^Jm  +  JhSlWySuJM)  = 

-  F'^WyS^JM  -  [/pw(-l)  -  urYWJiiJM  (3.45) 


A  new  matrix,  Kdu-,  can  be  defined  as: 

Kdu  =  {Jj^K^WyKiJM  +  JajW^Jm  +  JjjSyWySyJAf)  (3.46) 


Equation  3.45  can  now  be  rewritten  as: 


=  -K-^[F^WySyJM  -  [Ipu{-1)  -  utYWJC.JmY'  (3.47) 


Reinserting  the  terms  of  F  yields: 
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=  -KduliSxX  +  +  H^v  +  Had  -  rfWyS^JM- 

[4«(-l)  -  utYW^K^JmY  (3.48) 

One  final  reorganization  yields  the  equation  shown  below. 


I*  =  -K2^[-f'^WySuJM  +  {H^v-VHadfWySuJM+ 

WuK^Jm  +  Sl^WySuJM)  -  u'^WuKiJm  +  x'^S^WySuJMf  (3.49) 

This  equation  can  now  be  separated  into  the  following  matrices. 


Kr  =  [-WyS^JM]  e  U^*rty)xim*nu)  (3  5q) 

=  [HjWySM  €  (3  5^) 

Kd  =  [Hj WySuJM]  e  ^(P*^^+^Mrn*nu)  (3  52) 

Ku  =  [IJWuKiJm  +  Sl,WySM  e  (3.53) 

Kt  =  [-W^KiJm]  e  (3.54) 

=  [S^WySuJM]  e  3fJ^^x(m*nn)  (3  55) 


Kdu  =  [JLK^WuK^Jm  +  JLWaJm  +  JlS^WySM  e  (3  gg) 

Using  these  matrices,  z*  can  be  defined  as: 


z*{k)  =  -Kj^[f^{k)Kr  +  v'^{k)K^  +  d^ik)Kd+ 

u^{k  -  1)K^  +  ul{k)KT  +  x^{k)K:,f  (3.57) 

One  final  substitution  from  Equation  3.6  is  required  to  reach  the  optimal  input 
vector  Au*,  the  first  term  of  which  will  be  implemented  as  the  system  input  for  the 
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current  time  step. 


Au*^Jmz*  (3.58) 

3.3  Constrained  MPC  Optimal  Solution 

This  cost  function  can  also  be  minimized  subject  to  linear  equality  constraints  m 
the  inputs,  input-variations  and  outputs.  A  synopsis  of  these  constraints  are  shown 
below. 


< 

u{k  -k  i|fc)  < 

(3.59) 

Auf'' 

< 

Au{k  +  i\k)  <  Auf^^\  i  =  0,  •  •  • 

,p-l  (3.60) 

-€  + 

< 

y(k  -f-  i  +  l\k)  <  y^^^  +  e 

(3.61) 

Au{k  +  j\k) 

= 

o 

II 

(3.62) 

e 

< 

0 

(3.63) 

Constraints  on  the  inputs  and  input-variations  are  treated  as  hard  constraints, 
meaning  that  they  are  never  violated.  These  constraints  would  be  equivalent  to  the 
physical  limits  of  a  given  actuator.  To  avoid  infeasible  optimization  problems,  the 
output  constraints  are  treated  at  soft  constraints,  meaning  that  they  may  be  violated, 
if  required  to  maintain  optimization  feasibility.  This  is  allowed  by  an  addition  to  the 
cost  function  called  a  slack  variable,  e.  The  slack  variable  is  normally  equal  to  zero. 
However,  when  a  given  output  crosses  it’s  constraint  limit,  e  becomes  equal  to  the 
magnitude  of  violation,  e  is  then  squared,  multiplied  by  a  weighting  factor,  and 
added  to  the  function  cost,  which  was  shown  in  Equation  3.5.  Thus,  while  violating 
constraints  greatly  increases  the  cost  of  a  given  input  sequence,  the  optimization 
problem  will  never  become  infeasible. 

Using  the  previously  defined  matrices,  we  can  now  write  the  cost  function  to  be 
optimized,  including  the  slack  variable  term: 
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J  =  +  z^KduZ  +  2[f'^Kr  +  +  u^{k  -  l)Ku  +  uIKt  +  x^{k)K^]z  (3.64) 

Likewise,  the  constraints  can  be  realized  as: 


Equation  3.64  can  now  be  minimized,  subject  to  Equation  3.65  and  the  constraint 
e  >  0.  This  optimization  can  be  solved  using  the  optimizer  of  choice.  An  excel¬ 
lent  implementation  of  this  constrained  optimal  control  can  be  found  in  a  controller 
created  by  Dr.  Alberto  Bemporad  and  associates  from  the  Automatic  Control  Lab¬ 
oratory  in  Zurich,  Switzerland.  This  controller  uses  the  Matlab  optimizer  routine, 
DANTZGMP.m  to  solve  the  quadratic  problem  previously  shown.  For  a  sHghtly 
more  detailed  treatment  of  the  constrained  equation  setup,  see  [1].  Because  of  the 
limited  time  and  computational  resources  available,  this  research  will  deal  exclusively 
with  unconstrained  MFC  control  theory. 


3.4  Blocking  and  Basis  Functions 

As  mentioned  previously,  Jm  is  a  matrix  used  to  impose  additional  constraints  on 
the  optimum  Au{k).  This  matrix  can  impose  either  a  blocking  function  or  a  basis 


51 


function  on  the  optimum  solution.  The  matrix  Jm  decreases  the  complexity  of  the 
optimization  problem  while  still  allowing,  in  most  cases,  enough  freedom  for  the 
optimizer  to  solve  the  problem.  Utilizing  the  matrix  Jm  can  significantly  reduce  the 
number  of  rows  in  the  matrices  shown  in  Equation  3.57,  thus  lowering  the  calculation 
load  for  every  optimizer  cycle. 

Blocking  Functions 

To  further  constrain  the  optimization  problem,  Jm  can  impose  blocking  requirements 
on  the  problem.  This  means  that  the  terms  in  I  act  individually  to  constrain  the 
terms  of  Au  at  certain  times  during  the  projection  horizon.  To  illustrate  this  concept, 
an  example  is  offered.  Let  the  matrix  Jm  be  defined  as  shown  below.  In  all  cases, 
the  variable  m  is  equal  to  the  number  of  columns  in  This  would  require  that 
w(0)  =  'u(l),  u{2)  =  u{3)  =  u{4)  and  ^(S)  =  u(6). 


1 

o 

o 

Amo 

0  0  0 

Aui 

0 

0/0 

Au2 

-02 

0  0  0 

dictates  that 

Am3 

0 

0  0  0 

A'U4 

0 

o 

o 

Ams 

0  0  0 

Aue 

0 

Figure  3-3  depicts  a  possible  input  and  input  variation  schedule  which  could  result 
from  this  example  Jm- 
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<1 


0  1  2  3  4  5  6  7 

Prediction  Time  (k) 


Figure  3-3:  1st  Blocking  function  example:  Inputs  and  Input  Increments 


A  second  example,  showing  another  variation  of  blocking  functions  is  included 
below.  Again,  the  Jm  representing  this  function  is  shown,  as  well  as  the  input 
variation  schedule  which  it  would  produce.  This  Jm  matrix  would  required  that 
Am(0)  =  A«(l),  Au{2)  =  Au{3)  —  Au(4)  and  Att(5)  =  Att(6). 


7  0  0 

O 

<1 

Zi 

7  0  0 

Aul 

O 

o 

<! 

Z2 

0  7  0 

dictates  that 

Auz 

0  7  0 

A'U4 

Z2 

O 

o 

Am5 

Zz 

1 - 

o 

o 

1 _ 

1 

I> 

f 

1 _ 

.  . 

Figure  3-4  depicts  a  possible  input  and  input  variation  schedule  which  could  result 
from  this  blocking  function. 
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5 


^3 


qL _ I _ I - 1 - 1 - 1 - i - 1 

0  1  2  3  4  5  6  7 


Prediction  Time  (k) 

Figure  3-4:  2nd  Blocking  function  example:  Inputs  and  Input  Increments 


Basis  Functions 

Another  approach  to  constraining  the  optimization  problem  is  to  use  the  matrix  Jm 
as  a  basis  function.  This  means  that  the  terms  in  I  act  in  concert  to  constrain  the 
terms  of  Aw  during  the  projection  horizon.  To  illustrate  this  concept,  an  example 
is  offered.  Let  the  matrix  Jm  be  defined  as  shovm  on  page  55.  In  this  case,  each 
column  is  a  separate  basis  function.  The  terms  of  z  determine  the  weighting  placed 
on  each  function  in  the  resulting  Aw,  which  is  also  showm.  For  simplicity’s  sake,  we 
will  assume  that  the  number  of  inputs,  nw,  is  one. 
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Jm  = 


1 

0 

0 

Auo 

^1 

1 

0 

0 

Aui 

^1 

1 

0.5 

0 

Au2 

Zi  +  0,5z2 

1 

1 

0.33 

dictates  that 

t> 

s 

Zi  ~h  Z2  -i~  0.332^3 

1 

1 

0.66 

Au4 

^1  +  ^2  +  0.662^3 

1 

1 

1 

Aus 

2i  +  22  +  23 

1 

1 

1 

1 

< 

_ 1 

2i  +  22  +  23 

Because  of  the  integration  of  the  matrix  Jm  into  the  cost  function,  as  previously 
shown,  Jm  reduces  the  size  of  the  problem  linearly  from  size  p  to  size  m,  without 
sacrificing  significant  accuracy.  The  results  of  a  comparison  between  the  performance 
of  a  system  using  the  Jm  matrix  shown  above  and  that  of  a  system  using  an  identity 
Jm  matrix  of  size  p  are  shown  below.  An  identity  Jm  matrix  represents  a  system 
with  no  imposed  blocking  or  basis  function  constraints. 


Example  Jm 

Identity  Jm 

Process  Time  (sec) 

0.0200 

0.0500 

Function  Cost  (-) 

0.0364 

0.0321 

Table  3.1:  Jm  performance  comparison 


Utihzing  the  example  Jm  yields  an  increase  in  the  cost  function  J  of  approximately 
13.4%.  However,  it  also  decreases  the  processing  time  necessary  to  complete  one 
computation  cycle  by  60%.  This  trade  between  minimizing  the  cost  and  lowering  the 
processing  speed  exists  for  all  Jm  matrices,  with  the  most  optimum  solution,  in  terms 
of  cost,  being  that  offered  by  the  identity  Jm  matrix.  The  optimum  solution,  in  terms 
of  speed,  is  dependant  upon  the  control  task  requirements. 
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3.5  MPC  Example 


To  further  clarify  the  process  used  to  calculate  an  input  command  using  MPC,  an  ex¬ 
ample  is  offered.  The  example  discrete  state  space  system  will  be  drawn  from  the  K-1 
linearized  model  at  the  point  of  maximum  dynamic  pressure,  approximately  80  sec¬ 
onds  into  a  nominal  flight.  The  matrices  and  necessary  trajectories  are  shown  below. 
As  previously  stated,  no  estimated  unmeasured  disturbances,  d{k),  were  considered 
for  this  research.  The  previously  discussed  example  Jm  will  be  used,  indicating  a 
prediction  horizon  of  7. 

As  shown  below,  there  are  four  states,  \9  9  u  tc],  where  u  and  w  are  the  body 
velocities  in  the  X  and  Z  inertial  directions  respectively,  as  defined  in  Chapter  2  and 
9  is  vehicle  pitch.  The  input  is  the  pitch  gimbal  command,  5p  and  the  disturbance  is 
the  wind  velocity  in  the  negative  Z  direction.  The  output  is  the  angle  of  attack  (a) 
of  the  body. 


0 

0.4504 

-0.0001 

0.0003 

3.1761 

-0.0003 

1 

0 

0 

0 

Bu  = 

0 

Bx,  — 

0 

0 

33.1546 

-0.0106 

0.0042 

-25.5243 

-0.0042 

0 

-75.5361 

0.0047 

-0.0192 

57.0848 

0.0192 

C  = 


0  1  -0.0003  0.0007 


-0.0007 


In  addition,  the  input  weighting  matrices,  Wu  and  H  a,  are  both  identity  matrices 
of  size  p,  while  the  output  weighting  matrix,  Wy,  is  two  times  an  identity  matrix  of 
size  p. 

Next  we  will  define  the  reference  trajectory,  f{k),  the  disturbance  trajectory,  v{k), 
the  nominal,  or  trim,  output  trajectory,  ^(fc),  the  nominal  input  trajectory,  the 

nominal  state,  xo{k),  the  actual  state,  x{k)  and  the  most  recent  input,  u{k  —  1).  f{k) 
represents  the  commanded  a  trajectory.  v(k)  represents  the  wind  disturbance  which 
is  projected  to  affect  the  plant  at  each  time  step.  Notice  that  this  vector  is  one  term 
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longer  than  all  other  trajectory  vectors,  due  to  the  prediction  formulation.  yo{k), 
uo(k)  and  xo{k)  represent  the  trim  output,  input  and  state  values  that  were  used 
during  the  linearization  process.  u{k  —  1)  represents  the  gimbal  command  calculated 
and  implemented  during  the  previous  optimization  processor  cycle. 


r{k) 


r  ■  n 

151 

r  -1 

-0.0222 

149 

0 

0 

-0.0210 

147 

0 

0 

-0.0199 

144 

0 

0 

-0.0187 

v{k)  — 

151 

yo{k)  = 

0 

Uo{k)  = 

0 

-0.0175 

158 

0 

0 

-0.0164 

155 

0 

0 

-0.0146 

0 

0 

L.  J 

157 

L.  J 

0 

-0.0229 

Xo{k)  = 

-0.4204 

x{k)  = 

-0.3188 

1120 

1126 

500.6 

549.8 

u{k  —  1)  = 


-0.01474 


Using  these  variables,  the  following  matrices  and  vectors  were  created,  using  the 
previously  defined  procedures. 


1 

1 

0 

0 

0 

0 

0 

1 - 

o 

1 

1 

1 

0 

0 

0 

0 

0 

1 

1 

1 

1 

0 

0 

0 

0 

1 

II 

1 

1 

1 

1 

0 

0 

0 

1 

1 

1 

1 

1 

1 

0 

0 

1 

1 

1 

1 

1 

1 

1 

0 

^ — 1 

_ 1 

1 

1 

1 

1 

1 

1 

1 
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1.0000 

-0.0628 

0.0000 

-0.0000 

0.0476 

-0.0628 

0.4517 

-0.0001 

0.0003 

3.2227 

0.4517 

-0.0543 

0.0000 

-0.0000 

3.0429 

-0.0543 

0.2056 

-0.0000 

0.0001 

5.1  = 

4.4760 

0.2056 

-0.0362 

0.0000 

-0.0000 

4.3124 

-0.0362 

0.0943 

-0.0000 

0.0001 

4.9643 

0.0943 

-0.0217 

0.0000 

-0.0000 

4.8533 

0.0476 

0 

0 

0 

0 

0 

0 

3.2227 

0.0476 

0 

0 

0 

0 

0 

3.0429 

3.2227 

0.0476 

0 

0 

0 

0 

4.4760 

3.0429 

3.2227 

0.0476 

0 

0 

0 

4.3124 

4.4760 

3.0429 

3.2227 

0.0476 

0 

0 

4.9643 

4.3124 

4.4760 

3.0429 

3.2227 

0.0476 

0 

4.8533 

4.9643 

4.3124 

4.4760 

3.0429 

3.2227 

0.0476 

0.0147 

-0.7000 

0 

0 

0 

0 

0 

0 

-0.3003 

0.0147 

-0.7000 

0 

0 

0 

0 

0 

0.0250 

-0.3003 

0.0147 

-0.7000 

0 

0 

0 

0 

-0.1360 

0.0250 

-0.3003 

0.0147 

-0.7000 

0 

0 

0 

0.0191 

-0.1360 

0.0250 

-0.3003 

0.0147 

-0.7000 

0 

0 

-0.0621 

0.0191 

-0.1360 

0.0250 

-0.3003 

0.0147 

-0.7000 

0 

0.0122 

-0.0621 

0.0191 

-0.1360 

0.0250 

-0.3003 

0.0147 

-0.7000 

*le—3 
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-0.0952 

0 

0 

-6.5407 

0 

0 

-12.6264 

-0.0476 

0 

-21.5783 

-3.3179 

-0.0314 

-30.2032 

-9.5835 

-2.1898 

-40.1318 

-17.1024 

-6.3575 

-49.8383 

-25.8908 

-13.5114 

Ku=  755.4407  279.3801  113.0192 


K-jp  — 


-0.0059 

-0.0097 

-0.0146 

-0.0234 

-0.0255 

-0.0355 

-0.0274 

-0.0349 

-1.0000 
-2.0000 
-3.0000  - 
-4.0000  - 
-5.0000  - 
-6.0000  - 
-7.0000  - 


)  -0.0010 
J  -0.0025 
)  -0.0026 
I  -0.0060 
)  -0.0067 
)  -0.0142 
I  -0.0116 
)  -0.0181 

0 

0 

-0.5000 
-1.5000  - 
-2.5000  - 
-3.5000  - 
-4.5000  - 


)  -0.0002 
)  -0.0010 
)  -0.0006 
)  -0.0023 

r  -0.0016 

l  -0.0055 
)  -0.0043 
L  -0.0095 

0 

0 

0 

-0.3300 

-0.9900 

-1.9900 

-2.9900 


13.6717 

3.6333 

1.4921 

1 

2.9842 

1.2462 

0.5398 

8.3081 

1.3831 

0.2333 

II 

1.2462 

0.5783 

0.2662 

-0.0020 

-0.0003 

-0.0001 

0.5398 

0.2662 

0.1304 

0.0059 

0.0010 

0.0002 

I- 

-1 

*le3 


These  structures  were  used  to  calculate  z*{k)  (Equation  3.57),  which  was  used  to 
calculate  Au*(k)  (Equation  3.58),  which  was  used  to  calculate  y(k)  (Equation  3.30) 
and  u*{k)  (Equation  3.7),  which  were  used  to  calculate  the  function  cost,  J  (Equation 
3.31).  These  vectors  are  shown  on  page  60. 
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0.0261 

-0.1315 

0.0261 

-0.0466 

0.0261 

0.0040 

-0.0376 

-0.0442 

Au*(k)  = 

-0.0109 

y{k)  = 

0.0070 

0.0217 

-0.0037 

-0.0315 

0.0037 

-0.0188 

0.0037 

-0.0262 

u*{k) 


0.0114 

0.0375 

0.0416 

0.0307 

0.0269 

0.0306 

0.0342 


J- 


0.0364 
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Chapter  4 


MPC  Architecture  Development 


4.1  Architecture  Candidates 

Many  possible  control  layout  scenarios  exist  for  inserting  an  MPC  controller  into 
the  outer  and  possibly  inner  loops  of  a  rocket  control  system.  In  order  to  decide 
which  MPC  control  architectures  should  be  evaluated  in  this  thesis,  a  detailed  hnear 
evaluation  was  conducted  on  five  possible  control  layout  candidates.  Based  on  the 
results  of  this  analysis,  two  of  the  architectures  were  chosen  for  further  study.  The 
five  candidate  architectures  are  shown  below  and  on  the  following  pages,  along  with 
short  descriptions. 


Figure  4-1;  Candidate  Architecture  1:  MPC  outer  loop  with  a  controlhng  inner  loop 


The  first  architecture  is  referred  to  as  aSAS  and  consists  of  an  MPC  controller 
issuing  angle  of  attack  (a)  commands  to  an  inner  loop  stability  augmentation  system 
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(SAS),  which  then  issues  commanded  gimbal  angle  commands  (Sc)  to  the  gimbal 
actuator.  The  actuator  dynamics  then  issue  the  actual  gimbal  command  ((5p)  to  the 
vehicle. 


Figure  4-2:  Candidate  Architecture  2:  MFC  outer  loop  with  0  controlling  inner  loop 

The  second  architecture,  referred  to  as  6SAS,  is  very  similar  to  the  first,  except 
that  the  MFC  controller  issues  pitch  (0)  commands  to  an  inner  loop  SAS.  This  SAS 
then  issues  Sc  commands  to  the  actuator  dynamics,  which  in  turn  issue  the  true  Sp 
commands  to  the  vehicle. 


Figure  4-3:  Candidate  Architecture  3:  MFC  outer  loop  with  acceleration  direction 
controlling  inner  loop 


The  third  architecture,  referred  to  as  a  SAS,  consists  of  an  outer  loop  MFC  con¬ 
troller  issuing  acceleration  direction  commands  to  an  inner  loop  acceleration  direction 
controller  SAS.  In  turn,  this  SAS  then  issues  Sc  commands  to  the  actuator  dynamics, 
which  then  issued  the  final  Sp  command  to  the  vehicle. 
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Figure  4-4:  Candidate  Architecture  4:  MFC  outer/inner  loop 


The  fourth  architecture,  referred  to  as  NoSAS,  was  unique  among  the  five  con¬ 
troller  layouts  in  that  it  was  the  only  option  in  which  the  MFC  controller  performed 
the  duties  of  both  the  outer  guidance  loop  and  the  inner  stability  augmentation  loop. 
The  Sc  commands  were  issued  directly  to  the  actuator  dynamics,  which  then  issued 
the  actual  Sp  commands  to  the  vehicle. 


Figure  4-5:  Candidate  Architecture  5:  MFC  outer  loop  with  5  controlling  inner  loop 


The  fifth  architecture,  referred  to  as  SSAS,  consisted  of  an  MFC  controller  issuing 
gimbal  angle  commands  to  an  inner  loop  S  SAS.  This  SAS  then  issued  Sc  commands 
to  the  actuator  dynamics,  which  issued  the  actual  Sp  commands  to  the  vehicle. 

The  performances  of  all  five  candidate  architectures  were  compared  to  the  per¬ 
formance  of  a  baseline  load-relief  controller  layout.  This  controller  was  designed  to 
represent  a  current  launch  vehicle  load  relief  scheme,  as  detailed  by  Moreno  [5].  A 
schematic  of  this  baseline  controller  layout  is  shown  in  Figure  4-6. 
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Figure  4-6:  Baseline  Architecture  Layout 


The  baseline  load  relief  controller  was  designed  to  follow  a  nominal  pitch  trajec¬ 
tory.  If  unexpected  wind  disturbances  introduced  additional  a,  the  load  relief  channel 
would  inject  a  separate  signal  into  the  feedback,  which  would  act  to  bring  a  back  to 
zero.  Because  all  prospective  architectures  were  rated  relative  to  this  baseline,  the 
actual  performance  of  the  baseline  was  less  important  than  the  performance  of  the 
alternative  MFC  controllers  in  relation  to  the  baseline. 

The  controller  gains  were  chosen  to  produce  controllers  with  bandwidths  of  below 
2  rad/sec.  All  gains  are  unitless  except  Kd,  which  has  units  of  seconds.  The  gains 
used  in  these  comparisons  are  shown  below. 


aSAS 

esAS 

~aSAS 

NoSAS 

SSAS 

Baseline 

Ka 

1 

1 

0.3 

- 

1 

- 

Kd 

0.9 

0.9 

0.35 

- 

0.9 

0.9 

Ks 

0.4 

- 

Kp 

1 

1 

- 

- 

1 

0.7 

Kst 

1 

- 

- 

3 

- 

- 

Lr 

- 

- 

-  ! 

- 

- 

0.5 

Table  4.1:  Architecture  Gain  Values 


Because  all  vehicle  pitch  dynamics  were  represented  by  linear  state-space  models, 
analysis  was  conducted  only  around  the  flight  point  represented  by  the  linearized 
model.  This  point  was  chosen  to  be  the  point  of  maximum  dynamic  pressure,  which 
was  approximately  80  seconds  into  a  nominal  flight.  Because  the  vehicle  is  most 
aerodynamically  unstable  at  this  flight  point,  it  is  normally  considered  the  most 
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critical  design  point  in  a  launch  vehicle’s  ascent.  Simulations  were  restricted  to  a  ten 
second  increment  of  the  flight,  centered  at  the  design  point. 


4.2  Linear  Vehicle  Pitch  Dynamics 

The  linear  state-space  pitch  dynamics  model  used  for  this  analysis  was  produced  by 
the  same  process  used  to  create  the  plant  ID  information  at  every  flight  point  for  the 
non-linear  MFC  controller.  The  pitch  dynamics  represented  the  vehicle’s  behavior  at 
approximately  80  seconds  into  a  nominal  flight.  This  LTI  model  was  in  the  standard 
state  space  form  shown  in  Equation  3.1.  This  model  was  then  reduced  to  only  contain 
the  a,  $  and  0  states.  In  addition,  the  effects  of  a  wind  disturbance  on  the  states  and 
output  were  included  in  the  linear  model  via  the  appropriate  and  measured 
disturbance  matrices. 


4.3  Wind  Models 

Both  synthetic  and  actual  wind  samples  were  used  for  this  comparison.  Actual  wind 
profiles  were  sampled  at  the  appropriate  altitude  for  the  time  period  during  launch 
associated  with  the  vicinity  of  the  flight  point  to  produce  the  actual  wind  samples. 
Synthetic  wind  samples  were  produced  by  passing  a  random  signal  through  a  low- 
pass  filter  which  produced  a  simulated  wind  disturbance  profile.  This  wind  and  an 
appropriate  vehicle  speed  were  then  used  to  calculate  the  a  that  a  vehicle  flying 
through  this  wind  would  experience.  Every  random  profile  was  linked  to  a  unique 
random  number  or  “seed” ,  which  was  used  to  initialize  the  random  signal  generator. 


4.4  Architecture  Comparison  Results 

The  results  shown  in  Figure  4-7  represent  the  analysis  of  1000  synthetic  wind  profiles 
and  127  actual  wind  profiles.  Options  1-5  represent  the  aSAS,  OS  AS,  ~aSAS, 
NoSAS  and  SSAS  architecture  options,  respectively. 
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1  2  3  4  5 

Architecture  Candidate  Option  # 

Figure  4-7:  Architecture  Comparison  Results 


The  results  from  these  trials  indicate  that  the  6SAS  option  was  clearly  deficient 
in  both  a  reduction  and  bending  moment  reduction  as  compared  to  the  other  four 
options.  The  '^SAS  showed  acceptable  performance  in  a  reduction,  but  performed 
slightly  worse  than  the  remaining  three  options  in  bending  moment  reduction.  The 
NoSAS  option  showed  the  best  overall  reduction  of  bending  moments,  as  well  as 
excellent  a  reduction  capability.  Because  of  this  performance  and  its  unique  control 
layout,  this  architecture  was  selected  for  further  investigation.  In  addition  to  the 
non-SAS  approach,  a  inner-loop  SAS  approach  was  desired  to  serve  as  an  additional 
comparison  case.  While  the  remaining  two  architectures,  SSAS  and  aSAS,  showed 
nearly  identical  results,  the  aSAS  option  was  chosen  because  of  its  more  familiar 
and  intuitive  controller  layout.  Unlike  these  two  controller  architectures,  which  both 
featured  a  feedback,  the  9 SAS  and  ~a  SAS  options  seem  to  have  been  hampered  by 
the  lack  of  direct  a  feedback.  The  direct  control  of  a  was  considered  advantageous 
in  terms  of  both  controller  integration  and  trajectory  design  integration. 
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The  results  of  this  down-select  process  were  a  SAS  equipped  and  a  non-SAS 
equipped  MFC  control  layout.  These  two  architectures  would  allow  the  research 
to  investigate  the  advantages  and  disadvantages  of  each  type  of  system  in  relation 
to  the  current  load-relief  standard,  while  providing  further  insight  into  the  design 
process  that  applies  to  each  approach.  For  the  remainder  of  this  thesis,  the  MFC 
controller  with  an  inner-loop  a  controller  will  be  referred  to  as  AS  AS,  and  the  MFC 
controller  with  no  inner-loop  SAS  will  be  referred  to  as  NS  AS. 
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Chapter  5 


Wind  Discussion 

5.1  Wind  Profiles 

The  wind  profiles  used  for  this  research  were  gathered  from  the  Eastern  Test  Range 
(ETR),  located  at  Kennedy  Space  Center,  FL,  and  the  Western  Test  Range  (WTR), 
located  at  Vandenburg  AFB,  CA.  All  wind  profiles  were  gathered  using  jimsphere 
wind  observation  methods.  Reflective  aluminum  balloons,  such  as  the  one  shown 
below,  are  released  and  tracked  by  radar  as  they  ascend  through  the  atmosphere. 
The  radar  track  is  then  analyzed  to  calculate  the  wind  speed  at  each  altitude. 


Figure  5-1;  Jimsphere  Wind  Measurement  Balloon 
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The  ETR  wind  set  included  1,927  wind  profiles,  recorded  between  the  years  1964 
and  2000.  The  WTR  wind  set  included  1,437  wind  profiles,  recorded  between  1965  and 
1980.  Although  some  of  the  profiles  were  unusable,  due  to  tracking  and  transmission 
errors,  a  large  majority  of  the  wind  profiles,  1,141,  were  judged  to  be  accurate  and 
realistic  representations  of  the  actual  wind  state  at  the  time  of  the  observation.  The 
wind  data  gathered  from  both  ranges  included  wind  speed  and  wind  direction,  in 
roughly  80  foot  increments,  from  the  surface  to  an  altitude  of  approximately  95,000 
feet.  Because  all  research  was  conducted  in  the  pitch  plane  only,  it  was  necessary  to 
decompose  the  wind  velocity  into  North/South  and  East/West  components.  Because 
the  WTR  is  primarily  used  to  launch  payloads  into  polar  orbits,  all  tests  using  WTR 
wind  profiles  were  conducted  using  the  North/South  wind  component.  In  the  same 
way,  all  ETR  tests  were  completed  using  the  East/West  wind  component,  as  the  ETR 
is  primarily  used  to  launch  payloads  into  equatorial  orbits.  Figure  5-2  shows  a  single 
representative  wind  profile  from  the  ETR. 


Figure  5-2:  Sample  Wind  Profile 


Figure  5-3  shows  all  127  ETR  wind  profiles  used  in  the  MPC  controller  design 
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process.  This  set  will  hereafter  be  called  the  design  wind  set.  This  wind  set  is  a 
representative  sample  of  the  complete  wind  set  used  in  the  actual  system  comparisons. 

The  large  number  of  wind  profiles  capture  a  broad  cross-section  of  weather  types 
and  wind  behavior,  including  headwind  and  tailwind  conditions.  This  wide  range  of 
wind  conditions  presents  a  load  relief  problem  for  the  control  system.  Load  relief  is 
most  critical  during  portions  of  ascent  where  dynamic  pressure  (Q)  is  high.  During 
this  time,  all  MFC  load  relief  controllers  seek  to  minimize  vehicle  a  and  thus  Qa,  the 
product  of  Q  and  a.  Because  a  is  directly  affected  by  wind  speed  and  direction,  only 
one  unique  pitch  trajectory  exists  which  will  produce  the  minimum  a  for  any  given 
wind  profile  while  avoiding  excessive  Q  levels.  To  find  this  unique  pitch  trajectory 
for  each  wind  profile,  a  Matlab  Trajectory  Designer  was  created. 


Wind  Speed  (ft/sec) 

Figure  5-3:  ETR  Wind  Profile  Set 


5.2  Trajectory  Design 

Because  of  low  vehicle  speed  in  the  early  seconds  of  flight,  the  vehicle  must  be  guided 
by  a  pitch  controller  until  it  has  reached  a  safe  altitude  and  speeds  great  enough 
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to  require  a  control.  In  addition,  because  of  the  low  vehicle  speeds,  no  significant 
load  relief  benefits  are  gained  by  minimizing  a  at  low  speeds.  Because  of  this,  all 
launch  vehicle  trajectories  in  this  research  can  be  divided  into  two  major  portions: 
the  pitch  control  portion  and  the  a  control  portion.  To  quantify  each  trajectory 
further,  the  pitch  control  portion  can  be  divided  into  three  more  sections.  The  first 
section  consists  of  a  vertical  rise  from  the  launch  pad,  for  some  designated  rise  time, 
Tr-  After  this  initial  rise,  the  vehicle  enters  the  second  section,  in  which  it  execiites  a 
kick  maneuver  to  some  designated  pitch  angle,  kickO,  over  some  designated  kick  time, 
Tfe.  Following  this,  the  vehicle  enters  the  third  portion  of  the  flight,  where  it  flies  at 
the  constant  pitch  angle,  kickO,  until  a  given  transition  Q  level  is  reached.  At  this 
point,  the  vehicle  transitions  into  load-relief  control  and  enters  the  fourth  portion 
of  the  trajectory,  which  lasts  for  the  remainder  of  the  boost  phase.  Load  relief  is 
achieved  in  the  initial  trajectory  design  process  using  a  standard  PID  a  controller,  in 
the  form  illustrated  in  Figure  6-1. 

Each  acceptable  trajectory  must  meet  several  end  conditions.  The  K-1  has  a 
nominal  boost  phase  flight  time  of  134  seconds.  At  the  end  of  the  initial  boost 
phase,  vehicle  staging  requirements  dictate  that  the  vehicle  must  be  at  a  staging 
Q  of  approximately  60  Ibf/ft^  (psf)-  A  third,  implicit  requirement  states  that  the 
vehicle  must  not  experience  a.  Q  oi  more  than  600  psf  at  any  point  during  the  initial 
boost  phase.  The  Trajectory  Designer  works  to  meet  the  staging  dynamic  pressure 
constraint  by  iteratively  adjusting  the  parameters  Tk  and  kickO  until  a  134  second 
trajectory  is  found  which  produces  a  final  Q  that  falls  within  2  psf  of  the  staging 
requirement.  If  both  of  the  explicit  requirements  are  met,  the  implicit  maximum  Q 
limit  will  not  be  exceeded.  The  trajectory  design  process  is  illustrated  in  Figure  5-4. 

Because  the  wind  profiles  only  extend  to  altitudes  of  approximately  55,000  feet, 
each  profile  was  augmented  with  extra  wind  data  to  produce  a  smooth  transition  from 
the  final  wind  speed  to  zero,  as  seen  in  Figure  5-3.  However,  because  all  load  relief 
comparisons  are  completed  below  the  augmented  data,  this  contrived  wind  data  does 
not  present  a  problem  with  simulation  accuracy.  The  additional  wind  information  is 
used  only  by  the  Trajectory  Designer  and  has  no  effect  on  the  pitch  trajectory  during 
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Gravity  Turn 

a  =  o 

Load  Relief 
SEGMENT  4 


Constant  Attitude 
SEGMENT  3 


Kick  Maneuver 
SEGMENT  2 


Launch  Site 


Staging  Conditions 
Q  =  60  psf 
t  =  134  sec 


Trajectory  defined  in  terms  of  four  segments: 

•  SEGMENT  1:  Vertical  Rise 

•  SEGMENT  2:  Kick  Maneuver 

•  SEGMENT  3:  Constant  Attitude 

•  SEGMENT  4:  a  =  0:  Load  Relief 

Trajectory  Designer  seeks  to  meet  Staging 
Conditions  by  iterating: 
kAngle  (kicke) 

Kick  Duration  (T,^) 


Kick  Duration  =  T, 


Kick  Start  =  T, 


Vertical  Rise 
SEGMENT  1 


Figure  5-4:  Trajectory  Design  Framework 


the  load  relief  portion  of  the  flight  that  is  measured. 

Once  an  acceptable  trajectory  was  found,  the  pitch  profile  was  saved  for  later  use, 
as  well  as  the  appropriate  trajectory  parameters  and  other  important  profiles,  such 
as  Q.  Each  data  set  was  indexed  to  the  particular  wind  profile  to  which  it  applied. 
It  is  also  important  to  note  that,  while  the  trajectories  produced  by  the  Trajectory 
Designer  are  very  close  approximations  to  actual  KAC  pitch  trajectories,  they  are 
not  accurate  enough  to  be  used  for  actual  launch  operations.  However,  they  are 
qualitatively  identical  to  real  launch  trajectories  and  serve  as  valid  approximations 
for  this  research,  which  is  not  guidance  oriented,  but  rather  load  relief  oriented. 

Figure  5-5  shows  the  9  trajectory  from  a  sample  wind  profile.  Notice  the  distinct 
portions  of  the  flight,  marked  by  distinct  changes  in  pitch. 
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Figure  5-5:  Sample  Pitch  Trajectory 


Figure  5-6:  ETR  Pitch  Trajectory  Set 


Figure  5-6  shows  a  plot  of  all  the  pitch  trajectories  designed  for  the  design  wind 
set.  Notice  the  wide  range  of  Tk  and  kickO  variables  present  in  response  to  the  varying 
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wind  conditions  contained  within  the  wind  set. 

Figure  5-7  shows  a  plot  of  all  the  Q  profiles  resulting  from  the  pitch  trajectories 
shown  in  Figure  5-6.  Notice  the  convergence  of  the  profiles  to  60  psf  at  the  staging 
time  of  134  seconds. 


Figure  5-7:  ETR  Q  Profile  Set 


Equivalent  trajectories,  although  not  shown,  were  designed  for  the  remaining  1,800 
ETR  wind  profiles  and  1,437  WTR  wind  profiles. 

5.3  Time  Lapse  Discussion 

The  current  method  of  pre-launch  wind  measurement  involves  the  release  and  tracking 
of  multiple  jimsphere  balloons,  as  discussed  previously.  Delays  caused  by  slow  balloon 
ascension  and  data  processing  axe  significant.  The  standard  wind  profiles  that  are 
loaded  into  a  current  launch  vehicle  load-relief  system  were  usually  measured  at  least 
90  minutes  before  the  vehicle  is  scheduled  to  launch.  This  delay  introduces  large 
amounts  of  uncertainty  into  the  launch  decision  process,  and  can  often  result  in 
cancellations  or  delays  due  to  quickly  changing  wind  conditions.  To  simulate  this 
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delay,  the  time  stamp  of  each  wind  profile  was  analyzed  to  find  the  wind  sets  with 
time  differences  of  90  minutes,  ±  15  minutes.  The  actual  wind  profiles  from  these 
sets  were  used  to  simulate  a  traditional  balloon-based  wind  measurement  system 
integrated  with  an  MFC  controller  or  a  modern  traditional  load-relief  system.  Below 
is  a  sample  wind  pair  selected  fi’om  the  ETR. 


Figure  5-8:  Sample  ETR  Wind  Pair 


Some  current  load-relief  systems  do  not  use  the  actual  wind  profile  during  flight, 
but  rather  a  polynomial  fit  of  each  wind  profile.  Depending  on  the  system,  they 
assume  the  wind  follows  a  6th  or  15th  order  polynomial  fit,  and  then  design  a  trajec¬ 
tory  to  minimize  bending  moments  or  a  for  that  modified  wind  profile.  Representing 
these  systems  required  a  special  trajectory  design  process.  Depending  on  the  system 
fidelity,  the  Trajectory  Designer  would  replace  the  actual  wind  profile  with  a  poly¬ 
nomial  fit  of  the  appropriate  order,  and  then  complete  the  trajectory  design  process 
as  described  previously.  Thus  each  wind  pair  would  contain  three  different  pitch 
trajectories,  each  corresponding  to  a  different  fidelity  of  wind  profile  representation. 

During  time  lapse  system  comparison,  the  more  recent  profile  in  each  wind  pair 
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would  be  loaded  into  the  simulation  as  the  actual  wind,  while  the  appropriate  pitch 
trajectory  designed  for  the  earher  wind  profile,  depending  on  the  control  system  being 
simulated,  would  be  loaded  into  the  vehicle  control  system.  Below  is  a  comparison 
of  a  sample  wind  profile  and  its  6th  and  15th  order  polynomial  fits.  The  complete 
results  of  these  comparisons  are  included  in  Chapter  8. 


Figure  5-9:  Comparison  of  Actual  Wind  Profile,  6th  Order  Fit  and  15th  Order  Fit 
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Chapter  6 


Baseline  Controller  Development 


6.1  PID  Controller  Description 

All  PID  controllers  used  in  the  vehicle  simulations  followed  a  standard  layout.  This 
type  of  controller  was  chosen  for  several  reasons.  Primarily,  this  controller  type  was 
chosen  because  it  was  somewhat  similar  to  the  controller  used  in  the  actual  K-1 
launch  vehicle.  The  position  of  the  four  gains  also  allows  the  controller  to  handle 
derivative  rate  feedback  and  address  steady  state  error,  while  also  allowing  a  high 
degree  of  flexibility  in  the  design  process.  This  layout  was  used  in  all  simulations  to 
control  vehicle  yaw  attitude.  In  addition,  this  controller  layout  was  used  to  control 
vehicle  pitch  attitude  in  the  baseline  vehicle  simulation.  Finally,  this  layout  served 
as  the  inner-loop,  a  controller  layout  for  the  SAS-equipped  MPC  vehicle  simulation. 


Figure  6-1:  Standard  PID  Controller  Layout 
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Four  gains,  Ka,  Kd,  Kp  and  Ks,  were  employed  with  a  lead  compensator  and  an 
integrator  as  shown  in  Figure  6-1. 

6.2  Controller  Gain  Selection  Method 

The  gain  selection  process  was  considered  essential  in  building  a  fair  comparison  struc¬ 
ture  between  the  baseline  simulation  and  the  MPC-equipped  simulations.  The  main 
method  used  to  ensure  fair  controller  comparisons  between  drastically  different  con¬ 
troller  architectures  was  to  match  the  closed  loop  band  widths  of  the  controller /vehicle 
systems.  A  system  bandwidth  of  1.75  rad/sec  was  chosen  as  a  target  for  the  controller 
gain  selection  process.  Restricting  the  bandwidth  allows  higher  order  effects  such  as 
fuel  slosh,  bending  modes,  Tail-Wags-Dog  effects,  computational  delays  and  transport 
delays  to  be  neglected  without  compromising  the  validity  of  the  simulation.  Because 
the  vehicle  dynamics  change  significantly  throughout  a  standard  flight,  two  of  the 
four  controller  gains  must  be  dynamically  scheduled.  A  Matlab  script  file  was  cre¬ 
ated  to  calculate  the  gains  which  would  yield  the  appropriate  controller  bandwidth 
at  each  flight  point.  Flight  points  were  defined  at  one  second  increments  throughout 
the  flight. 

In  the  interest  of  process  speed,  a  Matlab  search  algorithm,  fminsearch.m,  was 
used  to  conduct  the  gain  selection.  While  holding  Ka  and  Ks  steady  at  1.4  and 
0.1  respectively,  the  algorithm  adjusted  the  values  of  Kd  and  Kp  within  a  dedicated 
SIMULINK  model  containing  the  controller  and  plant  pitch  dynamics.  During  each 
search  iteration,  this  file  performed  a  linearization  of  the  model  and  then  calculated 
the  phase  and  magnitude  response  of  the  combined  linear  model  to  a  range  of  input 
frequencies.  A  linear  interpolation  was  then  used  to  find  the  closed  loop  bandwidth, 
defined  as  the  frequency  at  which  the  closed  loop  magnitude  response  fell  below  -3 
dB.  This  process  was  repeated  imtil  the  error  between  the  system  bandwidth  and  the 
target  bandwidth  fell  below  an  acceptable  tolerance  level,  in  this  case  0.01  rad/sec. 
A  plot  of  the  resulting  gain  schedules  is  shown  in  Figure  6-2,  followed  by  a  plot  of 
the  resulting  closed  loop  system  bandwidth  at  each  flight  point,  in  Figure  6-3. 
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Baseline  Controller  CL  Bandwidt 


Flight  Point  (sec) 


Figure  6-3:  Closed  Loop  Bandwidth  of  Baseline  Controller  at  each  Flight  Point 


Figure  6-4;  Closed  Loop  Bode  Plot  of  Baseline  Controller  at  each  Flight  Point 

A  bode  plot  of  the  system  at  each  flight  point  is  shown  in  Figure  6-4.  Notice 
the  1.75  rad/sec  frequency,  where  all  magnitude  responses  converge  at  the  -3  dB 
line.  Although  not  explicitly  included  in  the  search  algorithm,  gain  and  phase  margin 
requirements  of  6  dB  and  30  deg,  respectively,  were  implicitly  imposed  on  the  resulting 
systems.  The  gain  and  phase  margin  of  the  closed-loop  system  is  shown  in  Figures 
6-5  and  6-6  respectively.  With  an  average  gain  margin  of  more  than  22  dB  and  an 
average  phase  margin  of  more  than  100  deg,  the  controller  far  exceeded  the  stability 
requirements  at  all  flight  points. 

While  no  explicit  system  performance  requirements  were  imposed  upon  the  gain 
selection  process,  the  performance  of  the  controller  at  each  flight  point  was  evaluated 
to  ensure  that  large  performance  gaps  did  not  exist  between  the  baseline  controller’s 
performance  and  standard  “acceptable”  performance  criteria.  Figure  6-7  shows  a 
plot  of  the  step  response  of  the  controller  at  each  flight  point.  The  average  percent 
overshoot  was  less  than  20%,  while  the  average  rise  time  was  about  1.3  seconds,  both 
well  within  acceptable  limits. 
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Figure  6-5:  Gain  Margin  of  Baseline  Controller  at  each  Flight  Point 


Figure  6-6:  Phase  Margin  of  Baseline  Controller  at  each  Flight  Point 

This  method  of  determining  controller  gains  produced  a  control  system  with  a 
narrowly  constrained  bandwidth  throughout  the  entire  flight  envelope.  This  allowed 
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Figure  6-7:  Step  Response  of  Baseline  Controller  at  each  Flight  Point 

fair  performance  comparisons  between  this  controller  and  other  controllers  that  were 
similarly  constrained,  while  still  maintaining  adequate  performance  and  minimum 
required  levels  of  stability. 

Because  the  yaw  controller  was  not  required  to  handle  any  disturbances  or  non-zero 
trajectory  commands,  and  no  analysis  was  conducted  in  the  yaw  plane,  fixed  controller 
gains  were  considered  sufficient  for  all  simulation  yaw  controller  gains.  These  gains 
were  set  equal  to  the  scheduled  pitch  controller  gain  values  at  the  flight  point  nearest 
maximum  Q. 
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Chapter  7 


MPC  Controller  Development 

7.1  Controller  Design 

The  design  process  for  a  MPC  controller  requires  defining  multiple  design  parameters. 
These  parameters  include  the  controller  cycle  rate,  the  prediction  step  size,  the  Wy, 
Wui  and  Wa  weighting  matrices,  the  Jm  matrix  and  the  prediction  horizon.  Presently, 
there  is  no  set  procedure  which  can  be  followed  to  design  an  MPC  controller,  as  there 
is  for  an  LQR  controller,  for  example.  This  chapter  will  attempt  to  define  such  a 
process,  by  describing  the  steps  which  were  used  to  set  each  of  the  parameters  in  the 
research  controllers.  While  still  requiring  engineering  intuition  and  experience,  this 
process  will  serve  to  build  a  basic  framework  or  road  map,  which  can  be  followed  when 
designing  any  MPC  controller.  In  addition,  the  actual  implementation  of  the  MPC 
controller  in  SIMULINK  will  be  described.  All  simulations  included  in  this  chapter 
were  conducted  using  the  design  wind  set,  as  detailed  in  Chapter  5. 

7.1.1  MPC  Controller  Range 

It  is  important  to  note  that  the  MPC  controller  was  only  designed  to  operate  between 
40  and  100  seconds  in  a  nominal  launch  profile.  As  described  in  the  Trajectory 
Design  section  of  Chapter  2,  the  first  portion  of  the  boost  profile  is  executed  using  a 
standard  pitch  controller,  following  a  pre-designed  pitch  profile.  At  a  certain  Q  level. 
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Figure  7-1:  Standard  Dynamic  Pressure  Profile 

control  transitions  from  the  standard  PID  pitch  controller,  described  in  Chapter  6, 
to  the  MPC  controller,  which  seeks  to  minimize  a.  Because  of  the  Q  levels  chosen 
as  transition  points,  this  control  hand-off  always  occurs  after  40  seconds  into  the 
flight.  To  save  time,  no  MPC  control  information  was  calculated  before  this  point. 
In  addition,  all  wind  profiles  used  as  disturbance  models  terminate  at  altitudes  of 
approximately  55,000  ft,  which  occurs  at  about  95  seconds  into  a  nominal  flight. 
At  this  altitude,  Q  has  fallen  below  the  critical  level,  due  to  a  rapid  decrease  in 
atmospheric  pressure.  Because  of  this,  terminating  the  trajectory  comparisons  at 
95  seconds  into  each  flight  was  not  considered  problematic.  Thus,  no  MPC  control 
information  was  calculated  after  100  seconds.  Figure  7-1  shows  a  plot  of  a  curve 
from  a  standard  boost  trajectory.  All  system  comparisons  were  completed  within  the 
high  Q  region  from  approximately  50  seconds  to  95  seconds. 
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7.1.2  MPC  Design  Process  Overview 

The  steps  used  to  design  the  AS  AS  and  NS  AS  MPC  controllers  are  shown  in  the 
following  list: 

1.  Choose  desired  system  bandwidth 

2.  Choose  magnitude  of  look-ahead  time 

3.  Choose  relative  gains  of  Wy,  Wu,  and  weighting  matrices 

4.  Choose  prediction  step  size 

5.  Choose  Jm  matrix  form 

6.  Choose  controller  rate 

This  design  process  will  not  exhaustively  cover  the  entire  design  space.  This 
approach  is  not  feasible,  due  to  the  extremely  high  number  of  design  parameters,  and 
the  wide  range  of  possible  values  for  each  parameter.  Instead,  this  process  will  focus 
on  a  sequence  of  steps  which  reduce  the  range  of  possible  controller  designs  in  an 
orderly  fashion,  while  providing  satisfaction  that  pertinent  design  possibilities  have 
not  been  overlooked. 

7.1.3  System  Bandwidth 

As  detailed  in  Chapter  6,  a  constant  system  bandwidth  of  1.75  rad/sec  was  chosen 
to  ensure  a  fair  comparison  between  the  MPC  controllers  and  the  baseline  PID  con¬ 
troller.  This  target  bandwidth  is  representative  of  the  control  system  bandwidths 
found  in  multiple  contemporary  launch  vehicles  and  allows  higher-order  dynamics  to 
be  neglected  without  compromising  the  validity  of  the  simulation.  The  bandwidth 
choice  affects  all  following  performance  metrics  and  controller  parameters,  and  may 
require  iteration  to  address  the  unique  characteristics  of  certain  launch  vehicles.  A 
discussion  of  the  effects  of  the  closed-loop  bandwidth  on  this  vehicle  is  included  at 
the  end  of  this  section. 
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7.1.4  Look-Ahead  Time  Magnitude 

The  next  step  of  the  design  process  is  to  choose  the  look  ahead  time  that  the  MFC 
controller  will  use  in  its  computations.  For  this  step,  the  input  weighting  matrices, 
Wu,  and  Wa,  were  set  to  identity  matrices,  leaving  only  the  output  matrix,  Wy,  to 
be  manipulated.  This  allowed  the  simple  manipulation  of  the  ratio  between  the  input 
and  output  weighting  matrices  to  affect  the  system  closed  loop  bandwidth. 

A  Matlab  search  algorithm,  fminsearch.m,  was  used  to  conduct  the  Wy  gain 
selection.  While  holding  the  W-^  and  H'a  gains  constant  at  unity,  the  algorithm 
adjusted  the  value  of  the  Wy  gain  within  a  dedicated  SIMULINK  model  containing  the 
controller  and  plant  pitch  dynamics.  During  each  search  iteration,  this  file  performed 
a  linearization  of  the  model  and  then  calculated  the  phase  and  magnitude  response  of 
the  combined  linear  model  to  a  range  of  input  firequencies.  A  linear  interpolation  was 
then  used  to  find  the  bandwidth,  defined  as  the  frequency  at  which  the  magnitude 
response  fell  below  -3  dB.  This  process  was  repeated  at  each  appropriate  flight  point 
until  the  error  between  the  system  bandwidth  and  the  target  bandwidth  fell  below 
an  acceptable  tolerance  level,  in  this  case  0.01  rad/sec.  This  process  was  completed 
for  both  MFC  control  architectures  over  a  range  of  prediction  step  sizes  and  look 
ahead  times.  Figures  7-2  to  7-7  show  plots  of  the  Wy  gain  schedules  for  three  selected 
prediction  step  sizes  for  each  MFC  architecture.  Look  ahead  times  of  1  to  10  seconds 
were  analyzed.  Figure  7-8  shows  all  the  resulting  system  bandwidths.  All  bandwidths 
fell  within  the  prescribed  tolerances. 

These  graphs  show  that  the  gain  schedules  are  very  different  for  look  ahead  peri¬ 
ods  of  1  to  3  seconds  than  they  are  for  look  ahead  periods  of  4  seconds  and  greater. 
Above  4  seconds,  the  gain  schedules  show  very  little  change  as  the  look  ahead  period 
is  increased.  This  convergence  behavior  implies  that  increasing  the  look  ahead  pe¬ 
riod  beyond  4  seconds  will  provide  very  little  change  in  performance.  However,  the 
computational  load  will  continue  to  rise  in  proportion  to  the  square  of  the  look  ahead 
period,  due  to  the  structure  of  the  MFC  matrices.  To  show  that  increasing  the  look 
ahead  period  beyond  4  seconds  does  not  yield  a  significant  benefit  in  system 
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Flight  Point  (sec) 


Figure  7-3:  Wy  Gain  Schedules  for  AS  AS,  Prediction  Step  =  0.6  sec 
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Flight  Point  (sec) 

Figure  7-6:  Wy  Gain  Schedules  for  NSAS,  Prediction  Step  =  0.2  sec 
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Figure  7-7:  Wy  Gain  Schedules  for  NS  AS,  Prediction  Step  =  0.5  sec 


1.8 


1.79 


O  1,77 
(U 


^1.76 

§  1.74 

CQ 

h4 

U  1.73 


1.7 


40 


60  70  80 

Flight  Point  (sec) 


100 


Figure  7-8:  Closed  Loop  Bandwidths  for  all  Controller/Look  Ahead  Combinations 


performance,  a  set  of  comparison  runs  was  completed  using  the  design  wind  set.  The 
controller  in  the  first  run  used  a  look  ahead  period  of  4  seconds,  while  the  second 
used  a  look  ahead  period  of  10  seconds.  For  this  trial,  the  prediction  step  size  and 
controller  rate  were  both  set  to  0.2  seconds,  and  the  Jmi  H  u  ^•iid  H'a  matrices  were 
set  to  identities.  The  results  for  this  trail  are  shown  in  Table  7.1. 


MPC  Controller  Type 

No  SAS 

No  SAS 

a  SAS 

a  SAS 

Look  Ahead  Period  (sec) 

4 

10 

4 

10 

Average  Max  a  (deg) 

0.36053 

0.35713 

0.35585 

Average  Max  Qq  (psf-deg) 

153.28 

151.31 

150.64 

Average  Max  Bending  Moment  (ft-lb) 

155687 

155628 

216323 

216139 

Average  Max  Gimbal  Angle  (deg) 

0.37005 

0.37026 

0.56549 

0.57146 

Table  7.1:  Performance  Comparison  of  MPC  controllers  with  varying  look  ahead 
periods 


More  than  doubling  the  prediction  horizon  from  4  seconds  to  10  seconds  has 
a  negligible  effect  on  the  metrics  shown  above,  usually  resulting  in  a  less  than  one 
percent  increase  or  decrease  in  a  given  variable.  The  fact  that  there  were  no  significant 
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differences  between  the  two  trials  shows  that  increasing  the  look  ahead  period  beyond 
4  seconds  does  not  yield  any  appreciable  increase  in  performance.  Because  of  this, 
the  look  ahead  period  for  both  controllers  was  set  to  4  seconds. 

7.1.5  Wy,  Ww)  and  Wa  Matrix  Weighting 

When  choosing  the  weighting  matrices  Wy,  Wu,  and  Wa,  it  is  important  to  recognize 
that  the  significant  design  factors  are  the  ratios  between  the  three  matrices,  and 
not  the  actual  values  of  each  matrix.  This  is  because  MFC  seeks  to  minimize  the 
magnitude  of  a  cost  function.  Thus,  one  matrix  may  be  set  to  unity,  leaving  only  two 
matrix  gains  to  be  manipulated  in  the  design  process,  without  sacrificing  any  design 
fidelity.  The  Wy  gain  will  continue  to  be  manipulated  as  described  previously  to 
match  the  system  closed-loop  bandwidth  with  a  given  target  bandwidth.  To  examine 
the  impacts  of  possible  Wu  and  Wa  gain  ratios,  a  subset  of  all  the  possible  gain  ratios 
deemed  to  span  a  significant  portion  of  the  design  space  was  examined.  For  each 
architecture,  the  Wu  gain  was  set  to  unity,  and  four  different  gains  were  applied  to 
the  Wa  matrix:  0,  1/2,  1  and  2.  Simulations  using  the  design  wind  set  were  conducted 
with  controllers  using  these  four  gain  ratios  to  discover  the  sensitivity  of  controllers 
with  each  gain  ratio  to  noise  in  the  wind  disturbance  measurement  process.  Figures 
7-9  to  7-13  contain  graphs  showing  the  sensitivities  of  each  architecture/gain  ratio 
combination  to  random  wind  disturbance  measurement  inaccuracies  with  increasing 
standard  deviations,  ranging  from  0  to  20  ft/sec.  A  wind  disturbance  measurement 
inaccmacy  was  defined  as  the  difference  between  the  wind  profile  used  to  calculate 
the  MFC  control  and  the  wind  profile  that  actually  affected  the  vehicle  during  flight. 
For  all  test  runs,  prediction  step  size  and  controller  rate  were  set  to  0.2  seconds  and 
look-ahead  time  was  set  to  the  previously  chosen  value  of  4  seconds.  In  addition,  the 
effect  of  changing  the  weighting  on  errors  as  a  function  of  their  location  within  the 
prediction  horizon  was  not  investigated.  All  errors,  regardless  of  their  distance  from 
the  present  time  step,  were  weighted  equally.  Thus,  each  weighting  matrix  contained 
only  one  non-zero  value,  which  was  repeated  at  all  points  on  the  matrix  diagonal. 

Several  interesting  trends  exist  in  the  data.  Both  architectures  exhibit  similar 
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figure  7-9:  Mean  Peak  a  Sensitivity  to  H  a  Weighting 


Figure  7-10:  Mean  Peak  Qq  Sensitivity  to  H^a  Weightinj 


x10' 


Figure  7-11;  Mean  Peak  Bending  Moment  Sensitivity  to  Wa  Weighting 


Figure  7-12:  Mean  Peak  Gimbal  Angle  Sensitivity  to  Wa  Weighting 
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Figure  7-13:  Mean  Total  Gimbal  Angle  Travel  Sensitivity  to  Weighting 


a  and  Qa  behavior  at  low  levels  of  wind  uncertainty.  However,  for  these  metrics, 
the  NS  AS  controller  shows  a  much  higher  sensitivity  to  increasing  amounts  of  wind 
uncertainty  than  the  AS  AS  controller,  as  shown  in  Figures  7-9  and  7-10.  This  may 
partially  be  due  to  the  fact  that  the  AS  AS  controller  is  using  its  internal  SAS  to  com¬ 
pensate  for  the  actual  winds,  while  the  NS  AS  controller  has  no  equivalent  immediate 
feedback.  Both  architectures  show  low  sensitivity  to  changes  in  the  W/\  weighting 
value  for  these  two  parameters.  This  trend  is  reversed  in  gimbal  activity,  as  shown  in 
Figures  7-12  and  7-13.  Both  controllers  increase  their  peak  gimbal  angle  by  roughly 
the  same  percentage,  as  wind  uncertainty  increases,  although  the  NS  AS  controller’s 
gimbal  response  is  generally  2/3  the  size  of  the  AS  AS  controller’s  response.  In  addi¬ 
tion,  the  zero  weighting  in  the  AS  AS  controller  seems  to  be  significantly  more 
sensitive  to  increases  in  wind  uncertainty  than  the  equivalent  weighting  in  the  NS  AS 
controller.  This  is  most  likely  due  to  the  internal  SAS,  whose  performance  is  directly 
impacted  by  the  uncertainty  of  the  wind  disturbance  signal,  via  the  a  feedback.  As 
can  be  seen  in  Figure  7-11,  the  higher  gimbal  angles  of  the  AS  AS  controller  more  than 
offset  the  bending  moment  reductions  created  by  lower  a  performance  to  produce  to- 
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tal  bending  moments  that  are  approximately  1/3  higher  than  the  equivalent  NS  AS 
bending  moments,  at  low  wind  uncertainty  levels.  At  higher  wind  uncertainty  levels 
the  a  and  Qa  sensitivity  of  the  NS  AS  controller  begins  to  bring  the  two  controllers 
to  parity  in  moment  performance.  Notice  the  excessive  bending  moments  caused  by 
the  zero  weighting  in  the  AS  AS  controller  at  higher  wind  uncertainties. 

Figures  7-14  and  7-15  show  the  individual  bending  moment  plots  for  each  con¬ 
troller  architecture,  which  were  used  to  more  clearly  show  the  response  of  each  con¬ 
troller. 


Figure  7-14:  Mean  Peak  Bending  Moment  Sensitivity  of  ASAS  Controller  to  Wa 
Weighting 

Based  on  these  results  and  the  stated  purpose  of  this  research,  which  was  to 
lower  bending  moments,  the  NS  AS  Wa  weighting  was  set  to  0  and  the  ASAS  Wa 
weighting  was  set  to  1/2.  Recent  studies  have  shown  that  current  LIDAR  wind 
measurement  technology  is  capable  of  reUably  supplying  wind  disturbance  information 
with  accuracies  of  at  least  4  ft /sec  or  better.  Thus,  controller  performance  in  this 
region  was  of  primary  importance.  However,  because  of  system  robustness  concerns, 
the  controller  performance  above  4  ft/sec  was  also  considered.  The  NS  AS  controller 
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Wind  Noise  Standard  Deviation  (ft/sec) 

Figure  7-15:  Mean  Peak  Bending  Moment  Sensitivity  of  NSAS  Controller  to  Ba 
Weighting 

showed  very  little  sensitivity  to  changes  in  the  B  a  weighting  over  the  entire  wind 
uncertainty  range.  The  weighting  value  of  0  showed  the  best  performance  below  wind 
uncertainties  of  4  ft/sec,  and  very  good  performance  relative  to  the  other  weighting 
values  as  uncertainty  increased  above  that  point.  The  AS  AS  weighting  of  0  was 
deemed  unacceptable  because  of  the  poor  performance  at  higher  wind  uncertainties, 
and  the  marginal  advantage  it  showed  over  other  non-zero  weighting  values  at  low 
wind  uncertainties.  The  weighting  value  of  1/2  exhibited  nearly  the  same  performance 
at  wind  uncertainties  below  wind  noise  standard  deviations  of  4  ft/ sec,  and  very  good 
relative  performance  beyond  that  point. 

It  is  important  to  note  that  the  chosen  ratios  are  only  the  best  of  the  ratios  which 
were  evaluated.  It  is  quite  possible  that  more  nearly  optimal  weighting  values  could 
be  found  if  more  ratios  were  investigated.  However,  for  this  research,  this  search 
was  considered  sufficient  to  uncover  the  significant  performance  trends  caused  by  Wa 
ratio  manipulation. 

These  plots  show  that  below  the  wind  uncertainty  standard  deviation  level  of 
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4  ft/sec,  wind  uncertainty  has  virtually  no  impact  on  a  MPC  controller’s  bending 
moment  performance.  As  stated  previously,  this  level  of  wind  sensor  accuracy  is  very 
realistic,  given  the  current  state  of  LIDAR  sensor  technology  development.  However, 
this  type  of  analysis  should  be  retained  in  the  controller  development  process,  as  other 
systems  will  react  differently  to  disturbance  measurement  uncertainty,  and  utilize 
sensors  with  varying  levels  of  fidehty. 

7.1.6  Prediction  Step  Size 

The  prediction  step  size  is  the  time  step  size  used  to  propagate  the  system  forward 
in  time.  If  this  variable  is  set  to  0.2  seconds,  and  the  prediction  horizon  is  set  to 
10,  the  controller  wiU  propagate  the  system’s  performance  2  seconds  into  the  future, 
in  0.2  second  increments.  Because  the  controller  operates  on  discrete  state-space 
models,  this  increment  must  be  used  to  convert  all  continuous  system  models  to 
discrete  system  models.  Shannon’s  Sampling  Theorem  dictates  that  the  frequency  of 
samphng  must  be  at  least  twice  the  frequency  of  the  fastest  unstable  pole  frequency 
of  the  system  being  controlled  [7]. 

Figure  7-16  shows  a  plot  of  the  fastest  unstable  pole  in  each  of  the  two  MPC 
architectures  being  considered,  at  each  flight  point.  The  unstable  high-frequency 
aerodynamic  pole,  which  is  visible  in  the  NSAS  architecture,  is  stabilized  by  the 
inner-loop  SAS  in  the  AS  AS  architecture.  Thus  the  highest  frequency  unstable  pole 
in  the  AS  AS  controller  is  the  low-frequency  gravity-induced  pole. 

Because  the  inner-loop  SAS  is  controlling  a,  and  not  an  inertially  based  parameter, 
such  as  pitch  angle,  this  low  frequency  unstable  pole  cannot  be  stabihzed.  This  pole 
represents  the  velocity  vector,  which  slowly  falls  as  the  vehicle  travels  along  a  gravity 
turn  trajectory.  Because  the  controller  is  attempting  to  follow  this  velocity  vector 
with  the  vehicle’s  body  x-axis  vector,  the  unstable  pole  results.  Since  the  ASAS 
architecture  has  a  much  lower  frequency  unstable  pole,  larger  prediction  step  sizes 
are  possible  than  with  the  NSAS  architecture.  Although  this  low  frequency  pole  is 
present  in  the  NSAS  architecture,  it  is  overshadowed  by  the  unstable  aerodynamic 
pole,  as  shown  in  Figure  7-16. 
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Using  Shannon’s  Sampling  Theorem  to  calculate  the  maximum  sampling  periods 
possible  with  these  two  control  architectures  produces  the  curves  shown  in  Figure  7- 
17.  These  curves  show  a  theoretical  maximum  prediction  step  size  of  approximately 
0.5  seconds  for  the  NSAS  architecture,  and  approximately  4  seconds  for  the  AS  AS 
architecture. 


Figure  7-16;  Comparison  of  Highest  Frequency  Unstable  Plant  Poles 

The  available  computational  power  and  controller  implementation  must  also  be 
considered  when  picking  the  prediction  step  size.  If  processor  speed  is  limited,  or 
the  optimizer  will  be  operating  in  real  time,  a  larger  step  size  may  be  necessary  to 
achieve  the  desired  look-ahead  distance  with  the  available  computational  resources. 
If  more  processor  power  is  available,  or  the  controller  optimization  computations  will 
be  completed  before  flight,  a  smaller  prediction  step  may  be  possible.  Decreasing 
sampling  time  will  also  allow  more  accurate  disturbance  sampling,  but  will  decrease 
the  amount  of  look-ahead  time  that  any  given  prediction  horizon  provides.  Likewise, 
increasing  the  sampling  period  will  provide  a  longer  look-ahead  time  for  any  given 
prediction  horizon  and  equivalent  controller  cycle  computation  load.  For  flexibility 
and  implementation  reasons,  possible  sampling  periods  were  limited  to  those  that 
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Figure  7-17:  Comparison  of  Maximum  Allowable  Sampling  Periods 

could  be  combined  into  one  second  increments,  such  as  0.1  seconds,  0.2  seconds,  0.5 
seconds,  etc. 

Due  to  computational  limitations,  all  MPC  matrices  used  in  this  research  were 
calculated  before  laimch,  avoiding  many  of  the  computational  issues  that  would  exist 
for  a  controller  seeking  to  solve  the  MPC  optimization  problem  in  real  time.  However, 
an  effort  was  still  made  to  keep  the  MPC  matrix  sizes  as  small  as  possible  while  main¬ 
taining  all  significant  system  performance.  To  test  the  sensitivity  of  each  controller 
architecture  to  changes  in  prediction  step  size,  simulation  runs  were  accomplished 
using  the  design  wind  set  and  a  range  of  prediction  step  sizes.  Each  architecture  was 
tested  at  four  step  sizes.  The  ASAS  controller  was  tested  at  0.1,  0.2,  0.5  and  1.0 
second  prediction  step  sizes,  while  the  NSAS  controller  was  tested  at  0.05,  0.1,  0.2, 
and  0.5  second  step  sizes.  Step  sizes  of  greater  than  1  second  were  deemed  too  large 
to  accurately  sample  the  disturbance  input,  while  step  sizes  of  less  than  0.05  seconds 
were  deemed  too  small,  as  they  would  require  prediction  horizons  of  greater  than  80 
to  meet  the  previously  chosen  look-ahead  period  requirement  of  4  seconds.  Previously 
chosen  values  for  the  Wa  weighting  values  were  used.  The  results  of  these  tests  are 
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Figure  7-20:  Mean  Peak  Bending  Moment  Sensitivity  to  Prediction  Step  Size 
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Figure  7-21:  Mean  Peak  Gimbal  Angle  Sensitivity  to  Prediction  Step  Size 
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Figure  7-22:  Mean  Total  Gimbal  Angle  Travel  Sensitivity  to  Prediction  Step  Size 


As  in  previous  analysis,  the  a  and  Qq  curves  are  proportionally  similar,  with  the 
NS  AS  controller  exhibiting  a  bowl  shaped  response  as  the  prediction  step  size  is 
decreased.  The  increases  beyond  step  sizes  of  0.2  seconds  are  most  likely  due  to  a 
decrease  in  the  gimbal  activity  at  the  same  step  size  points.  The  AS  AS  controller 
follows  a  generally  decreasing  trend  for  both  of  these  metrics  as  step  size  is  decreased, 
with  a  slight  increase  at  the  smallest  step  size,  for  the  same  reasons  as  discussed 
above.  The  moment  curves  for  both  controllers  show  a  general  downward  trend  for 
the  evaluated  step  sizes.  However,  as  step  size  continues  to  decrease,  both  show 
progressively  smaller  reductions  in  bending  moments,  with  a  significant  levelling  of 
the  NS  AS  curve  below  step  sizes  of  0,2  seconds.  The  AS  AS  controller  shows  similar, 
though  less  drastic,  levelling  behavior.  The  gimbal  activity  plots  show  the  same 
general  trend  of  decreasing  gimbal  deflection  in  both  controllers  as  the  step  size  is 
reduced.  Gimbal  travel  follows  this  general  trend,  although  the  AS  AS  controller 
exhibits  an  interesting  increase  in  gimbal  travel  at  a  step  size  of  0.5  seconds. 

Based  on  these  results,  the  stated  goal  of  lowering  bending  moments,  and  the  desire 
to  avoid  excessive  computational  loads,  the  prediction  step  size  for  both  architectures 
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was  set  to  0.2  seconds.  In  each  case,  the  marginal  decreases  in  bending  moments 
beyond  this  point,  as  compared  to  earlier  reductions,  were  not  considered  significant 
enough  to  warrant  increasing  the  prediction  horizon  with  a  further  doubUng  from  20 
to  40. 

7.1.7  Jm  Matrix 

As  previously  stated,  all  MPC  matrices  used  for  this  research  were  calculated  before 
flight,  due  to  processor  speed  constraints.  Because  of  this,  the  Jm  matrix  used  for  all 
research  runs  was  set  equal  to  an  identity  matrix  of  the  appropriate  size.  In  general,  if 
a  non-identity  Jm  matrix  is  desired  for  processor  speed  reasons,  an  evaluation  process 
similar  to  those  described  in  previous  sections  of  this  chapter  should  be  conducted. 
The  performance  of  multiple  Jm  matrices  would  be  tested  until  one  was  found  that 
yielded  suitable  time  savings  while  still  meeting  a  given  set  of  controller  accuracy 
requirements. 

7.1.8  Controller  Rate 

Once  the  prediction  step  size  has  been  selected,  the  controller  cycle  rate  can  be  chosen. 
The  controller  rate  is  the  number  of  times  per  second  that  the  MPC  controller  cycles. 
If  this  variable  is  set  to  10  Hz,  the  controller  will  calculate  a  new  input  for  the  plant 
10  times  per  second.  Because  the  MPC  controller  implements  the  first  term  of  the 
optimal  input  sequence,  the  controller  rate  must  be  greater  than  or  equal  to  one 
over  the  prediction  step  size.  A  controller  rate  below  this  limit  would  result  in  the 
controller  implementing  a  control  input  for  longer  than  it  was  designed  to  be  used, 
which  would  compromise  the  stability  of  the  controller.  Increases  in  this  parameter 
may  be  beneficial,  but  are  limited  by  the  available  computing  power.  To  determine 
if  any  benefits  could  be  gained  from  increasing  the  controller  rate  size  beyond  the 
standard  rate,  a  series  of  trials  were  completed  using  the  design  wind  set.  Each 
controller  used  the  settings  previously  discussed,  with  controller  rates  of  5,  10  and  20 
Hz  for  both  architectures.  Previously  chosen  values  for  Wa  a.nd  prediction  step  size 
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Figure  7-25:  Mean  Peak  Bending  Moment  Sensitivity  to  Controller  Rate 
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Figure  7-26:  Mean  Peak  Gimbal  Angle  Sensitivity  to  Controller  Rate 
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Figure  7-27:  Mean  Total  Gimbal  Angle  Travel  Sensitivity  to  Controller  Rate 


These  plots  show  the  MFC  controller’s  relatively  low  level  of  sensitivity  to  in¬ 
creases  in  the  controller  rate.  Both  controllers  show  only  marginal  decreases  in  bend¬ 
ing  moments,  gimbal  deflection  and  gimbal  travel  above  rates  of  10  Hz.  Slight  in¬ 
creases  in  the  NS  AS  a  and  Qa  are  evident  as  controller  rate  increases.  The  AS  AS 
controller  shows  almost  no  response  to  increasing  controller  rate  in  the  a  and  Qa 
metrics. 

Based  on  these  results  and  the  stated  purpose  of  this  research,  which  was  to 
lower  bending  moments,  the  NS  AS  and  AS  AS  controller  rates  were  both  set  to  10 
Hz.  Neither  architecture  showed  benefits  above  this  rate  significant  enough  to  justify 
doubling  the  rate  from  10  to  20  Hz. 

7.1.9  Bandwidth  Manipulation  Sensitivity 

When  compared  to  the  industry  standard  baseline  controller,  the  two  MFC  controllers 
presented  in  this  chapter  show  excellent  load  reduction  characteristics,  as  detailed  in 
Chapter  8.  However,  these  lower  a,  Qa  and  bending  moment  metrics  are  gained 
at  the  price  of  significantly  increased  peak  gimbal  angle  deflection  and  total  gimbal 
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angle  travel.  Some  launch  vehicles  power  their  engine  gimbals  with  energy  bleeds  from 
the  combustion  chamber,  and  can  handle  large  increases  in  gimbal  activity  without 
problem.  However,  the  engine  gimbal  systems  of  many  modern  launch  vehicles  are 
powered  with  open-loop  pressure  systems,  which  are  only  capable  of  rotating  the 
engines  through  a  preset  number  of  degrees  before  running  out  of  power.  For  these 
vehicles,  the  high  gimbal  activity  shown  by  the  previously  detailed  MFC  controllers 
could  be  problematic. 

The  abihty  to  trade  increased  peak  ct,  Qa  and  bending  moment  metrics  for  de¬ 
creased  gimbal  activity  would  be  very  useful  for  launch  vehicles  operating  under 
these  limitations.  This  can  be  achieved  by  increasing  the  weighting  value  in  the  input 
weighting  matrix,  Wu-  This  will  increase  the  cost  associated  with  any  given  gimbal 
angle  movement,  which  will  naturally  decrease  the  control  systems  ability  to  respond 
to  wind  disturbances,  resulting  in  increased  peak  a,  Qa  and  bending  moment  met¬ 
rics.  As  mentioned  earlier  in  this  section,  the  important  variable  is  not  the  actual 
value  of  each  weighting  matrix,  but  rather  the  ratio  between  each  set  of  weighting 
values.  Thus  an  increase  in  the  weighting  value  of  Wy,  is  equivalent  to  an  appropri¬ 
ate  decrease  in  the  weighting  value  on  the  output  matrix,  Wy.  Lowering  this  value 
causes  a  decrease  in  the  closed-loop  system  bandwidth.  To  investigate  the  trade-off 
between  gimbal  activity  and  the  three  metrics  of  a,  Qa  and  bending  moments,  AS  AS 
and  NS  AS  controllers  with  bandwidths  of  1,  1.25,  1.5  and  1.75  rad/sec,  respectively, 
were  created  using  the  previously  chosen  variable  settings.  The  performance  of  these 
controllers  in  response  to  the  design  wind  set  is  shown  in  Figures  7-28  to  7-32. 

The  results  reveal  several  insights  into  the  sensitivity  of  the  MFC  controllers.  The 
AS  AS  controller  shows  only  very  minor  changes  in  its  performance  metrics  as  the 
bandwidth  is  decreased  from  1.75  to  1  rad/sec.  If  the  launch  vehicle  can  handle  the 
increased  gimbal  activity  that  the  AS  AS  controller  produces,  this  stationary  behavior 
can  be  an  advantage,  since  the  system  bandwidth  can  be  lowered  substantially  without 
significantly  affecting  the  load  reduction  performance  of  the  controller.  This  relative 
insensitivity  to  changes  in  bandwidth  is  due  to  the  a  SAS,  which  operates  at  very  high 
cycle  rates,  resulting  in  higher  gimbal  activity  than  the  NS  AS  controller,  regardless 
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of  system  bandwidth. 


MPC  Controller  Bandwidth 

Figure  7-28:  Angie  of  Attack  Comparison  for  Varying  Controller  Bandwidths 
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Figure  7-29:  Qa  Comparison  for  Varving  Controller  Bandwidths 
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Figure  7-30:  Bending  Moment  Comparison  for  Varying  Controller  Bandwidths 
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Figure  7-31:  Gimbal  Angle  Deflection  Comparison  for  Varying  Controller  Bandwidths 
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MPC  Controller  Bandwidth 

Figure  7-32:  Total  Gimbal  Angie  Travel  Comparison  for  Varying  Controller  Band- 
widths 

However,  if  the  ability  to  trade  peak  bending  moment  or  Qa  performance  in 
exchange  for  decreased  gimbal  activity  is  desired,  the  NS  AS  controller  must  be  used. 
As  shown  in  the  previous  figures,  the  NS  AS  controller  metrics  vary  significantly  as 
the  system  bandwidth  is  lowered  from  1.75  to  1  rad/sec.  These  changes  result  in  an 
approximate  50%  decrease  in  the  gimbal  angle  deflection  and  total  gimbal  angle  travel, 
in  exchange  for  an  approximate  35-40%  increase  in  a,  Qa  and  bending  moments  over 
this  range.  This  is  a  very  important  capability  if  special  consideration  must  be  taken 
for  a  given  launch  vehicle’s  engine  gimbal  movement  limitations. 

The  K-1  engine  gimbals  are  powered  by  an  accumulator  drawing  power  directly 
from  the  vehicle  and  are  not  subject  to  movement  limitations  which  approach  the 
levels  shown  in  this  thesis.  Because  of  this,  the  K-1  can  handle  large  increases  in 
engine  gimbal  movement  without  any  significant  penalty  to  the  vehicle  itself.  Thus, 
pursuant  to  the  stated  goal  of  this  research,  w^hich  is  to  minimize  bending  moments 
during  ascent,  the  previously  chosen  closed-loop  system  bandwidth  of  1.75  rad/sec 
will  be  used  during  the  remainder  of  this  research. 
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7.1.10  Controller  Variable  Summary 

The  results  of  the  controller  tuning  process  are  summarized  in  the  Table  7.2: 


MPC  Controller  Type 

No  SAS 

a  SAS 

Closed-Loop  System  Bandwidth  (rad/sec) 

1.75 

1.75 

Look  Ahead  Period  (sec) 

4 

4 

Prediction  Horizon  (-) 

20 

20 

Wy  Weighting  (-) 

Adjusted  for 

Bandwidth 

Adjusted  for 

Bandwidth 

Weighting  (-) 

1 

1 

Wa  Weighting  (-) 

0 

0.5 

Prediction  Step  Size  (sec) 

0.2 

0.2 

Jm  Matrix 

Identity 

Identity 

Controller  Rate  (Hz) 

10 

10 

Table  7.2:  MFC  Controller  Parameter  Setting  Summary 


7.2  Wind  Processing  Analysis 

As  described  in  Chapter  5,  the  wind  profiles  used  for  this  research  contain  wind  infor¬ 
mation  in  approximately  80  ft  increments.  Current  LIDAR  wind  sensors  are  capable 
of  discerning  wind  information  in  increments  of  between  approximately  35  to  245 
ft.  From  this  information,  it  is  possible  to  pre-process  the  wind  data  in  a  number 
of  different  fashions,  including  higher  order  polynomial  fits,  varying  sampling  incre¬ 
ments  or  moving  average  smoothing.  In  this  case,  moving  average  smoothing  was  the 
method  of  choice  because  it  retains  a  relatively  high  amount  of  accuracy,  compared 
to  sampling,  as  the  averaging  distance  is  increased.  Sampling  at  various  increments 
and  linearly  interpolating  between  data  points  was  also  investigated  and  found  to  be 
deficient,  due  to  poor  accuracy  retention  as  the  increment  distance  was  increased. 
Polynomial  fits  were  unacceptable  because  they  could  not  accurately  represent  wind 
profiles  at  very  high  levels  of  fidelity.  To  determine  whether  an  advantage  might  be 
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gained  from  some  level  of  wind  processing,  as  well  as  to  investigate  each  controllers 
response  to  less  and  less  accurate  wind  information,  both  MPC  architectures  were 
flown  against  the  design  wind  set,  using  forward  and  backward  moving  average  dis¬ 
tances  of  0,  1,  2,  4,  8,  16,  24,  and  32  data  points,  respectively.  The  spacing  between 
these  points  ranges  from  0  ft  to  approximately  2600  ft  in  each  direction,  hereafter 
called  the  ±  distance. 

Figures  7-48  to  7-52  show  an  interesting  set  of  trades  between  gimbal  activity 
and  bending  moments  for  both  MPC  controllers.  At  certain  points  along  the  curves, 
favorable  trades  exist  which  could  be  used  to  modulate  the  controller’s  performance  to 
compensate  for  individual  launch  vehicle  limitations  in  gimbal  movement  capability. 
As  in  the  bandwidth  manipulation  example,  the  AS  AS  controller  does  not  show 
a  significant  gimbal  activity  decrease  in  response  to  changes  in  sensed  wind  profile 
smoothing,  relative  to  the  much  larger  increases  in  the  a  and  Qa  metrics.  This  results 
in  bending  moment  metrics  which  increase  at  unacceptable  rates  in  response  to  small 
or  negligible  decreases  in  gimbal  activity.  As  noted  before,  this  is  because  of  the 
internal  a  SAS,  which  operates  separately  at  nearly  real-time  rates,  in  response  to  the 
actual,  unprocessed  wind  profile.  Thus,  the  gimbal  activity  is  relatively  unresponsive 
to  smoothing  of  the  sensed  disturbance  vector,  as  this  vector  is  loaded  into  the  MPC 
controller,  but  not  the  a  SAS. 

The  NS  AS  controller  shows  more  gimbal  sensitivity  to  the  smoothing  of  the  dis¬ 
turbance  vector.  Although  the  bending  moment  curve  is  continually  increasing,  there 
are  points  where  the  percentage  increase  in  bending  moments  is  significantly  out¬ 
weighed  by  the  percentage  decrease  in  gimbal  angle  deflection  and  total  gimbal  angle 
travel.  In  a  vehicle  with  limited  gimbal  movement  capability,  this  option  can  prove 
important.  An  example  smoothing  average  ±  distance  of  330  ft  shows  an  approxi¬ 
mately  3%  increase  in  bending  moments  in  exchange  for  an  approximate  decrease  in 
gimbal  angle  deflection  and  total  gimbal  angle  travel  of  20%. 

In  light  of  these  results,  the  stated  purpose  of  this  research,  and  the  K-l’s  robust 
ability  to  absorb  increases  in  commanded  total  gimbal  angle  travel,  all  wind  distur¬ 
bance  profiles  will  be  loaded  into  the  MPC  controllers  at  the  highest  level  of  fidelity 
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possible,  which  is  in  approximately  80  ft  increments  in  this  case. 


Moving  Average!  Distance  (ft) 


Figure  7-33:  Angle  of  Attack  Comparison  for  Varying  Wind  Speed  Moving  Average 
±  Distances 


Moving  Average!  Distance  (ft) 


Figure  7-34:  Qa  Comparison  for  Varying  Wind  Speed  Moving  Average  ±  Distances 
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Moving  Average!  Distance  (ft) 

Figure  7-35:  Bending  Moment  Comparison  for  Varying  Wind  Speed  Moving  Average 
±  Distances 


Figure  7-36:  Gimbal  Angle  Deflection  Comparison  for  Varying  Wind  Speed  Moving 
Average  ±  Distances 
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Moving  Average ±  Distance  (ft) 

Figure  7-37:  Total  Gimbal  Angle  Travel  Comparison  for  Varying  Wind  Speed  Moving 
Average  ±  Distances 


7.3  Controller  Stability  Analysis 

The  stability  analysis  of  each  MFC  controller  was  complicated  by  the  presence  of  a  low 
frequency  unstable  pole,  as  noted  previously.  This  pole  was  created  by  commanding 
the  vehicle  to  follow  a  variable  trajectory  that  was  not  inertially  based,  in  this  case 
a  given  a  trajectory.  Commanding  the  vehicle  to  follow  a  variable  such  as  pitch 
would  eliminate  this  pole.  However,  because  of  the  load-relief  focus  of  this  research, 
it  was  necessary  to  command  an  a  trajectory,  as  opposed  to  a  pitch  or  flight-path 
angle  trajectory.  This  mode  diverged  very  slowly,  and  was  very  unresponsive  to  gain 
manipulation.  To  calculate  the  gain  and  phase  margin  of  the  system,  a  dedicated 
SIMULINK  model  was  used,  much  like  the  model  used  to  calculate  closed-loop  system 
bandwidth.  To  calculate  the  gain  margin,  the  input  command  from  the  controller  to 
the  plant  was  multiplied  by  larger  and  larger  gains,  until  one  of  the  system  poles, 
aside  from  the  small  pole  due  to  oc  following,  crossed  the  juo  axis.  This  value  was  then 
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converted  to  find  the  gain  margin  of  the  system  in  dB.  The  phase  shift  was  calculated 
in  the  same  manner.  The  input  command  was  affected  by  larger  and  larger  time 
delays  until  the  first  stable  pole  was  driven  into  the  right  hand  plane.  The  delay 
was  implemented  by  a  2nd  order  pade  linear  approximation  of  the  actual  non-linear 
delay,  This  time  delay  was  then  used  with  the  magnitude  of  the  pole  crossover 

frequency  to  calculate  the  phase  margin.  This  process  was  repeated  at  each  flight 
point  for  both  the  AS  AS  and  NS  AS  controllers.  A  picture  of  the  SIMULINK  model 
used  to  calculate  the  gain  and  phase  margins  for  the  NS  AS  controller  is  included 
below. 


Figure  7-38:  Stability  Analysis  Model  for  NSAS  Controller 


Although  not  explicitly  included  in  the  design  process,  gain  and  phase  margin  re¬ 
quirements  of  6  dB  and  30  deg,  respectively,  were  implicitly  imposed  on  the  resulting 
systems.  With  an  average  gain  margin  of  more  than  18  dB  and  an  average  phase 
margin  of  more  than  55  deg,  both  MFC  controllers  far  exceeded  the  stability  require¬ 
ments  at  all  flight  points.  Figures  7-39  and  7-40  show  plots  of  sample  pole  tracks  as 
gain  and  time  delay  are  increased,  respectively.  The  gain  and  phase  margins  of  both 
controllers  are  shown  in  Figures  7-41  and  7-42. 

A  bode  plot  of  each  system  at  each  flight  point  is  shown  in  Figures  7-43  and  7-44. 
Notice  the  1.75  rad/sec  frequency,  where  all  magnitude  responses  converge  at  the  -3 
dB  line.  A  DC  gain  of  less  than  one  is  visible  in  the  AS  AS  bode  plot.  This  will  result 
in  slight  tracking  errors,  which  must  be  compensated  for  during  flight.  The  solution 
to  this  problem  is  discussed  in  the  Controller  Implementation  section. 
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Imaginary  Axis  (rad/sec' 


Real  Axis  (rad/sec) 


Figure  7-39:  Pole  Positions  for  Increasing  Gain  Values 


Real  Axis  (rad/sec) 

Figure  7-40:  Pole  Positions  for  Increasing  Time  Delay  Values 


Figure  7-43:  Closed  Loop  Bode  Plot  of  AS  AS  Controller  at  each  Flight  Point 


Figure  7-45:  Step  Response  of  ASAS  Controller  at  each  Flight  Point 


Figure  7-46:  Step  Response  of  NS  AS  Controller  at  each  Flight  Point 

While  no  explicit  system  performance  requirements  were  imposed  during  the  de¬ 
sign  process,  the  performance  of  each  controller  at  each  flight  point  was  evaluated 
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to  ensure  that  large  performance  gaps  did  not  exist  between  each  controller’s  perfor¬ 
mance  and  standard  “acceptable”  performance  criteria.  Figures  7-45  and  7-46  show 
plots  of  the  step  response  of  each  controller,  respectively,  at  each  flight  point.  The 
AS  AS  controller  average  percent  overshoot  was  less  than  5%,  while  the  average  rise 
time  was  about  1.4  seconds.  The  NS  AS  controller  average  percent  overshoot  was  less 
than  5%,  while  the  average  rise  time  was  about  1.6  seconds.  All  metrics  were  well 
within  acceptable  limits.  The  steady-state  error  of  the  AS  AS  controller  is  clearly 
visible,  and  a  very  small  steady  state  error  is  visible  for  most  of  the  NS  AS  controller 
step  responses. 

7.4  Controller  Implementation 

Chapter  3  detailed  the  construction  of  the  MFC  matrices:  Kr,  K^,  Kd,  Ku,  Kd-,  Kt 
and  Kdu.  Using  these  matrices,  the  optimal  input  parameter  sequence,  z*{k),  can  be 
calculated.  Because  the  Jm  matrix  used  for  this  research  was  an  identity  matrix,  this 
vector  is  equal  to  the  optimal  input  control  change  sequence,  Au*{k).  This  equation 
is  shown  below: 


Au*{k)  =  -Kj^[f'^{k)Kr  +  v'^{k)K^  +  <f{k)Kd+ 

u^{k  -  l)Ku  +  Ux{k)KT  +  x^{k)KxY  (7-1) 


This  equation  can  be  rewritten  as: 


Au{k)  =  -K^^[T'^{k)KrY'  -  K^^[v^{k)K,Y  -  K^^[d^{k)KdY- 

Kdu[^^{k  -  1)K^Y  -  K2::[ul{k)KTY  -  K-^[x^{k)KxY  (7-2) 

For  this  research,  only  measured  disturbances  were  considered,  thus  eliminating 
the  Kd  matrix.  The  remaining  five  terms  were  reorganized  to  create  the  following 


123 


matrices. 


K7-  = 

(7.3) 

Kv  = 

-KauKv 

(7.4) 

Ku  = 

(7.5) 

Kt  = 

(7.6) 

Kx  = 

(7.7) 

Using  these  newly  defined  terms,  Equation  7.2  can  be  rewritten  to  yield  the  opti¬ 
mal  control  sequence,  as  shown  below. 

Au  (k)  =  Kr  *  r(/c)  -f-  Kv  *  v(k)  -\-  Ku  *  u(k  —  1)  -I-  Kt  *  Mr(A^)  -I-  Ax  *  x(k)  (7.8) 

Because  only  the  first  term  of  any  given  optimal  control  sequence  calculated  by 
the  MFC  controller  is  implemented,  only  the  first  row  of  each  of  these  matrices  must 
be  stored  at  each  flight  point.  During  flight,  a  linear  interpolation  is  completed  be¬ 
tween  the  appropriate  flight  points  to  determine  the  correct  interpolated  K  vector. 
Once  each  of  these  vectors  is  calculated,  only  a  simple  series  of  vector  multiplica¬ 
tions  and  summations  is  required  to  calculate  the  first  term  of  the  optimal  control 
sequence.  This  allows  extremely  fast  execution  times  wdth  very  low  computational 
power  requirements.  Figure  7-47  shows  a  schematic  of  the  closed-form  SIMULINK 
MFC  controller  setup  used  by  both  MFC  architectures. 

The  summing  junction  immediately  before  the  controller  output  adds  every  suc¬ 
cessive  Au  term  to  the  last  commanded  input,  creating  the  next  commanded  input. 
To  implement  this  in  a  SIMULINK  environment,  the  unit  delay  is  required.  Each 
of  the  five  K  vectors  are  linearly  interpolated  in  the  blocks  shovm.  The  disturbance 
trajectory  is  a  vector  of  the  anticipated  wind  disturbance  magnitudes  at  each  predic¬ 
tion  step  from  the  present  time  to  the  prediction  horizon.  These  vectors  are  gathered 
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Unit  Delay 


Figure  7-47:  MFC  Controller  SIMULINK  Implementation 


from  a  look-up  table  relating  altitude  to  wind  disturbance  magnitude.  Because  this 
analysis  is  completed  only  in  the  pitch  plane,  the  wind  magnitude  information  is  in 
the  inertial  X-Z  plane.  The  anticipated  altitude  of  the  vehicle  is  calculated  using  a 
simple  kinematics  equation  and  the  vehicle’s  present  position,  velocity  and  accelera¬ 
tion.  The  state  vector  is  the  current  state  from  the  vehicle,  in  this  case  the  current 
[9  9  u  w]^  where  u  and  w  are  the  body  velocities  in  the  X  and  Z  inertial  directions 
respectively,  as  defined  in  Chapter  2.  YO,  UO  and  XO  represent  the  trim  output,  input 
and  state  values  that  were  present  during  the  linearization  process.  XO  is  the  current 
trim  state,  while  YO  and  UO  are  projections  of  the  trim  values  forward  in  time  to  the 
projection  horizon.  The  vector  aT  represents  the  commanded  output,  in  this  case 
some  commanded  a  trajectory.  All  sampling  times  in  the  input  ports,  output  ports 
and  digital  clocks  are  set  to  the  controller  rate,  while  the  prediction  step  is  used  in  all 
blocks  which  are  used  to  find  predicted  values  from  look-up  tables,  such  as  the  trim 
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states  and  a  trajectory. 


7.5  Controller  Bias  Processing  Analysis 

Due  to  the  MFC  implementation  method,  both  controllers  exhibited  tracking  errors 
during  flight.  This  caused  the  controller  to  drift  slightly  during  flight,  usually  resulting- 
in  an  a  trajectory  which  contained  a  1/2  to  1/4  degree  error,  in  relation  to  the 
commanded  zero  a  trajectory.  The  solution  to  this  problem  was  to  run  the  simulation 
twice.  The  first  run  would  serve  as  an  initialization  run,  while  the  second  run  would  be 
considered  the  actual  flight.  During  the  load  relief  portion  of  the  initialization  run,  the 
commanded  a  trajectory  would  be  set  to  zero,  and  the  resulting  a  trajectory  would 
be  recorded.  During  the  second  simulation  run,  the  commanded  a  trajectory  would 
be  set  equal  to  some  filtered  version  of  the  negative  of  the  previously  recorded  actual 
a  trajectory.  This  would  effectively  cancel  out  all  significant  drift  that  the  controller 
may  have  exhibited  during  the  initialization  run.  This  correction/command  signal  is 
represented  by  the  aT  vector  shown  in  Figure  7-47. 

Because  the  recorded  a  signal  contains  significant  high  frequency  content,  it  may 
be  advantageous  to  process  this  signal  before  using  it  as  the  commanded  a  trajectory. 
As  previously  stated  in  the  wind  processing  analysis  section,  the  processing  method 
of  choice  is  to  smooth  the  alpha  bias  signal  using  a  moving  average.  To  determine 
how  both  MFC  controllers  respond  to  different  amounts  of  smoothing,  both  MFC 
simulations  were  flown  against  the  design  wind  set  with  moving  average  ±  distances 
of  0,  2,  4,  8,  16,  24,  32,  40,  48,  56  and  64  data  points,  respectively.  These  equate  to 
dz  time  spans  ranging  from  0  to  6.4  seconds.  The  results  are  shown  in  Figures  7-48 
to  7-52. 

These  figures  show  an  interesting  set  of  metric  curves  from  the  MFC  controllers 
in  response  to  different  amounts  of  a  bias  signal  smoothing.  In  this  case,  it  seems 
possible  to  conduct  an  acceptable  trade  between  gimbal  activity  and  bending  moments 
for  both  controllers.  The  AS  AS  controller  shows  a  significant  response  to  smoothing 
of  this  signal,  as  it  is  a  primary  driver  of  the  gimbal  activity  produced  by  the  a  SAS. 
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■  0  1  2  3  4  5  6 

Moving  Average!  Time  Span  (sec) 

Figure  7-48;  Angle  of  Attack  Comparison  for  Varying  Bias  Moving  Average  ±  Time 
Spans 


0  1  2  3  4  5  6 

Moving  Average+  Time  Span  (sec) 

Figure  7-49:  Qa  Comparison  for  Varying  Bias  Moving  Average  ±  Time  Spans 
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Moving  Average!  Time  Span  (sec) 


Figure  7-50:  Bending  Moment  Comparison  for  Varying  Bias  Moving  Average  ±  Time 
Spans 
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Figure  7-51:  Gimbal  Angle  Deflection  Comparison  for  Varying  Bias  Moving  Average 
±  Time  Spans 
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Figure  7-52:  Total  Gimbal  Angle  Travel  Comparison  for  Varying  Bias  Moving  Average 
±  Time  Spans 


The  NS  AS  controller  also  shows  a  very  strong  trade  relationship  between  gimbal 
activity  and  bending  moment  performance.  One  example  of  this  trade  is  at  the 
moving  average  ±  time  span  of  4  seconds.  At  this  point,  an  increase  of  just  over  1% 
in  bending  moments  produces  an  approximate  decrease  of  26%  in  total  gimbal  angle 
travel  and  30%  in  gimbal  angle  deflection.  These  types  of  gimbal  activity  reductions 
are  quite  significant  and  could  be  very  useful  in  the  pursuit  of  lowered  gimbal  activity 
for  launch  vehicles  operating  with  MFC  control  systems. 

This  chapter  has  examined  three  separate  methods  of  trading  bending  moment 
performance  for  gimbal  activity.  These  methods,  controller  closed-loop  bandwidth 
manipulation,  wind  disturbance  signal  processing  and  bias  correction  signal  process¬ 
ing,  can  be  implemented  individually  or  in  combination  with  each  other,  depending 
on  the  vehicle  restrictions  present  in  each  situation.  Although  each  method  has  been 
explored  individually,  no  combination  studies  have  been  conducted,  which  may  yield 
new  insights  into  the  problem  of  high  input  activity  for  MFC-controlled  systems. 


129 


[This  page  intentionally  left  blank. 


130 


Chapter  8 


Simulation  Comparison  Results 


A  wide  range  of  comparison  trial  runs  were  completed  to  fully  evaluate  the  perfor¬ 
mance  of  both  MPC  controllers  and  the  baseline  controller.  As  described  in  Chapters 
1  and  5,  the  space  laimch  industry  currently  relies  upon  a  balloon-based  system  to 
measure  the  winds  during  launch  operations.  This  balloon  measurement  process  yields 
wind  information  that  is  an  average  of  90  minutes  old  at  the  time  of  launch.  During 
pre-launch  preparations,  this  wind  information  is  checked  against  a  set  of  go/no  go 
criteria  to  ensure  that  certain  wind  hmits  are  not  violated.  In  addition,  several  mod¬ 
ern  load  rehef  systems  use  a  polynomial  representation  of  the  wind  in  this  pre-launch 
process  to  design  a  trajectory  which  minimizes  certain  vehicle  metrics,  such  as  ol  and 
Qa.  The  vehicle  then  follows  this  trajectory  during  flight,  using  a  standard  pitch 
controller.  The  earliest  forms  of  these  load  relief  systems  used  3rd  order  fits  to  model 
the  wind.  More  recent  upgrades  to  these  systems  utilize  6th  and  15th  order  fits. 

8.1  Full  Wind  Set  Comparisons 

A  comprehensive  wind  set  of  over  3,000  wind  profiles  was  used  for  the  controller 
comparisons.  To  simulate  a  time  lapse  between  wind  measurement  and  the  actual 
laimch,  the  wind  sets  were  analyzed  and  pairs  of  wind  profiles  occurring  90  minutes 
apart,  ±  15  minutes,  were  selected.  38  of  these  pairs  were  found  in  the  ETR  wind  set, 
and  156  pairs  were  found  in  the  WTR  wind  set.  During  simulation,  the  vehicle  would 
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follow  a  commanded  trajectory  designed  in  response  to  the  earlier  wind  of  the  pair, 
while  the  later  wind  would  be  applied  to  the  vehicle  during  flight.  All  comparisons 
involving  time  delays  were  conducted  using  these  wind  pairs. 

The  benefits  of  only  implementing  an  MFC  controller  or  a  LIDAR  wind  sensor 
individually  were  also  investigated.  To  evaluate  the  possible  benefits  of  using  a  LI¬ 
DAR  wind  sensor  alone  without  upgrading  the  on-board  vehicle  load-relief  system, 
simulations  were  conducted  against  all  wind  profiles  for  each  controller,  using  the 
trajectory  designed  to  minimize  a  for  that  particular  wind.  It  is  important  to  recall 
that  while  the  baseline  controller  follows  this  pre-planned  trajectory  for  the  full  dura¬ 
tion  of  the  flight,  the  MFC  controllers  transition  into  an  a  following  load-relief  mode 
upon  passing  a  Q  level  of  200  psf.  To  evaluate  the  possible  benefits  of  using  a  MFC 
controller  alone  without  improving  the  wind  measurement  equipment  or  process,  the 
MFC  simulations  were  flown  using  the  time  lapse  wind  pairs  described  previously. 
While  the  later  wind  profile  would  affect  the  vehicle  during  flight,  the  earlier  wind 
profile  would  be  loaded  into  the  MFC  controller’s  prediction  process. 

Finally,  to  show  the  benefits  of  an  MFC  controller,  utilizing  wind  information 
provided  by  a  LIDAR  wind  sensor,  both  MFC  controllers  were  flown  against  all  wind 
profiles  while  using  current  wind  information  for  the  controller’s  prediction  process. 
The  data  resulting  from  these  simulation  runs  is  summarized  in  the  following  pages. 
The  five  metrics  used  in  previous  sections,  peak  a,  peak  Qa,  peak  bending  moments, 
peak  gimbal  angle  defleetion  and  total  gimbal  angle  travel  are  shown,  for  the  wind  sets 
from  each  test  range.  Delay  and  no  delay  results  are  shown  for  five  different  controller 
and  wind  fidelity  combinations,  with  a  small  graphical  offset  between  each  for  added 
clarity.  Both  the  average  and  standard  deviation  of  each  metric  at  each  point  are 
shown.  The  five  combinations  include  the  baseline  controller  following  trajectories 
designed  using  the  actual  wind  profile,  denoted  as  Inf,  as  well  as  trajectories  designed 
using  15th  order  and  6th  order  wind  fits.  The  remaining  two  combinations  are  the 
AS  AS  and  NS  AS  MFC  controllers,  following  actual  wind  profiles,  both  present  and 
delayed. 
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Figure  8-1:  Controller  Peak  a  Comparison  for  ETR 
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Figure  8-2:  Controller  Peak  a  Comparison  for  WTR 
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Figure  8-3:  Controller  Peak  Qa  Comparison  for  ETR 


Figure  8-4:  Controller  Peak  Qa  Comparison  for  WTR 


Figure  8-5;  Controller  Peak  Bending  Moment  Comparison  for  ETR 


Figure  8-6:  Controller  Peak  Bending  Moment  Comparison  for  WTR 
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%  Reduction  using  AS  AS  Controller  with  Actual,  No  Delay  ETR  Winds: 

Base  6  D 

Base  15  D 

Base  Inf  D 

Base  6  ND 

Base  15  ND 

Base  Inf  ND 

a 

75.79 

74.33 

76.60 

72.95 

66.07 

67.82 

Qa 

75.94 

74.38 

76.84 

72.42 

65.13 

67.42 

B.M. 

71.45 

70.10 

74.62 

67.99 

59.95 

67.25 

G.A.D. 

-156.59 

-124.25 

15.98 

-195.35 

-234.55 

8.80 

T.G.A.T. 

-528.53 

-503.26 

-25.59 

-520.57 

-520.22 

-34.63 

%  Reduction  using  AS  AS  Controller  with  Actual,  No  Delay  WTR  Winds: 

Base  6  D 

Base  15  D 

Base  Inf  D 

Base  6  ND 

Base  15  ND 

Base  Inf  ND 

a 

73.01 

71.05 

74.18 

70.27 

62.87 

64.59 

Qa 

72.77 

71.10 

74.17 

69.17 

61.80 

64.49 

B.M. 

68.80 

66.92 

72.54 

65.37 

56.97 

64.82 

G.A.D. 

-170.49 

-175.77 

10.80 

-213.00 

-270.14 

5.90 

T.G.A.T. 

-568.35 

-561.82 

-33.06 

-561.30 

-567.85 

-33.79 

%  Reduction  using  NS  AS  Controller  with  Actual,  No  Delay  ETR  Winds: 

Base  6  D 

Base  15  D 

Base  Inf  D 

Base  6  ND 

Base  15  ND 

Base  Inf  ND 

a 

75.49 

74.02 

76.31 

72.62 

65.65 

67.42 

Qa 

75.82 

74.25 

76.72 

72.28 

64.95 

67.25 

B.M. 

78.72 

77.71 

81.08 

76.14 

70.15 

75.59 

G.A.D. 

-106.37 

-80.36 

32.42 

-137.55 

-169.08 

26.65 

T.G.A.T. 

-283.56 

-268.14 

23.36 

-278.71 

-278.49 

17.84 

%  Reduction  using  NS  AS  Controller  with  Actual,  No  Delay  WTR  Winds: 

Base  6  D 

Base  15  D 

Base  Inf  D 

Base  6  ND 

Base  15  ND 

Base  Inf  ND 

a 

72.79 

70.82 

73.97 

70.02 

62.57 

64.30 

Qa 

72.90 

71.24 

74.28 

69.31 

61.98 

64.65 

B.M. 

76.83 

75.44 

79.61 

74.29 

68.05 

73.88 

G.A.D. 

-115.94 

-120.15 

28.79 

-149.88 

-195.50 

24.88 

T.G.A.T. 

-305.12 

-301.16 

19.35 

-300.85 

-304.82 

18.90 

Table  8.1:  Metric  %  Reductions  Due  to  MFC  Controllers  and  LIDAR  Wind  Sensors 


for  Both  Test  Ranges 
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Table  8.1  shows  %  reductions  between  each  MFC  controller  and  the  baseline  con¬ 
troller  under  various  wind  measurement  and  baseline  controller  configuration  com¬ 
binations.  Because  a  percent  reduction  is  represented  by  a  positive  number,  a  neg¬ 
ative  number  represents  a  percent  increase  from  the  baseline  value.  The  complete 
results  from  these  trials  can  be  found  in  Appendix  C.  Notice  that  no  analysis  was 
conducted  using  the  MFC  controllers  in  combination  with  processed  wind  profiles. 
Using  smoothed  wind  data,  as  described  in  Chapter  7,  can  produce  positive  results  in 
certain  situations,  depending  on  the  individual  launch  vehicle  and  its  gimbal  move¬ 
ment  capabilities.  However,  these  gains  cause  significant  increases  in  the  a,  Qa  and 
bending  moment  metrics.  For  the  K-1,  these  issues  do  not  exist.  Therefore,  in  order 
to  achieve  maximum  reduction  in  bending  moments,  actual  wind  data  and  bias  signals 
were  used  for  these  results.  A  more  productive  method  of  lowering  gimbal  activity 
without  sacrificing  a  significant  amount  of  bending  moment  performance  is  to  smooth 
the  bias  signal,  as  shown  at  the  end  of  Chapter  7.  A  set  of  bias  smoothing  results 
are  included  in  the  next  section.  Baseline  controller  configurations  are  represented  by 
6,  15  and  Inf,  which  are  the  order  of  the  polynomial  fits  used  to  represent  the  wind 
in  each  controller,  respectively.  Balloon  and  LIDAR  wind  measmement  methods  are 
represented  by  D  (delay)  and  ND  (no  delay),  respectively. 

Many  interesting  observations  can  be  made  after  studying  these  figures  and  tables. 
Significant  areas  of  interest  are  the  standard  deviations  of  each  controller /wind  fidelity 
combination  at  different  time  delays,  the  gimbal  behavior  of  the  MFC  controllers,  as 
compared  to  that  of  various  versions  of  the  baseline  controller,  and  the  impact  of  the 
LIDAR  wind  sensor  and  the  effect  of  wind  fidelity  on  all  controllers’  performances. 

The  first  fact  to  note  is  that  in  all  cases,  both  MFC  controllers  are  able  to  reduce 
the  bending  moments  sustained  by  the  launch  vehicle  during  ascent  by  greater  than 
50%,  ’when  mated  with  a  LIDAR  wind  sensor  which  supplies  current  wind  data  to 
the  controller.  This  achievement  is  made  even  more  significant  when  the  standard 
deviation  of  the  peak  bending  moment  metric  is  taken  into  account.  Both  MFC 
controllers  show  bending  moment  standard  deviations  which  are  approximately  four 
times  smaller  than  equivalent  delayed  baseline  results.  This  increase  in  performance 
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certainty  would  allow  the  vehicle  to  fly  safely  in  rapidly  changing  weather  conditions 
without  increased  risk,  since  the  less  conservative  launch  decision  process  would  be 
offset  by  increased  confldence  in  the  vehicle  performance  capability. 

The  benefits  of  the  LIDAR  wind  sensor  mated  with  any  of  the  simulation  con¬ 
trollers  are  obvious.  When  the  LIDAR  wind  sensor  replaces  the  balloon  measurement 
system  in  conjunction  with  a  baseline  controller,  measurement  uncertainty  is  reduced 
significantly.  This  is  reflected  in  the  decrease  in  a,  Qa  and  bending  moment  standard 
deviations  of  approximately  50%.  In  general,  the  more  advanced  the  control  system, 
the  more  sensitive  it  is  to  wind  measurement  errors.  This  is  especially  evident  in  the 
case  of  the  MFC  controllers,  which  perform  extremely  well  with  current  wind  infor¬ 
mation,  but  much  worse  with  delayed  wind  information.  Relative  to  their  nominal, 
no  delay,  performance,  the  MFC  controllers  actually  perform  worse  than  the  baseline 
controllers  with  delayed  wind  information.  In  contrast,  the  baseline  controller  using 
a  6th  order  fit  of  the  wind  is  affected  the  least  by  delayed  information,  since  the  low 
order  fit  acts  as  a  filter.  This  filtering  behavior  removes  all  the  high  frequency  wind 
information,  and  leaves  only  the  general  wind  speed  trends,  which  are  much  less  likely 
to  change  significantly  over  a  90  minute  period.  As  the  wind  profile  fidelity  increases, 
errors  between  the  actual  wind  and  the  sensed  wind  used  in  the  controller  become 
more  and  more  detrimental  to  system  performance. 

The  gimbal  behavior  of  each  controller  is  also  very  interesting.  The  baseline  con¬ 
trollers  using  6th  and  15th  order  wind  fits  show  extremely  low  gimbal  angle  deflections 
and  total  gimbal  angle  travel,  because  of  the  smooth  simulated  wind  profiles.  This 
allows  the  trajectory  designer  to  produce  commanded  pitch  trajectories  which  require 
very  little  gimbal  activity  to  follow.  However,  when  the  actual  wind  profile  is  used, 
the  pitch  profile  produced  by  the  trajectory  designer  is  very  jagged  with  high  fre¬ 
quency  characteristics.  When  the  baseline  pitch  controller  attempts  to  follows  this 
commanded  trajectory,  it  requires  much  more  gimbal  activity,  as  shown  in  Figures 
8-7  through  8-10.  However,  because  the  vehicle  is  using  a  pitch  controller  to  mini¬ 
mize  a,  the  extra  gimbal  activitj^  does  not  translate  efficiently  into  lower  a  and  Qa 
performance.  In  contrast,  the  MFC  controllers  are  able  to  translate  increased  gimbal 
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activity  into  lower  a,  Qa  and  bending  moment  performance,  since  they  operate  in 
a  direct  load-relief  fashion,  where  the  variable  being  controlled  is  directly  related  to 
aerodynamic  loads. 

Finally,  it  is  important  to  realize  that  by  design,  MFC  controllers  produce  better 
results  vs.  a  traditional  control  system  as  the  wind  profiles  become  more  and  more 
erratic  and  extreme.  Because  of  the  wide  range  of  wind  profiles  included  in  this  work, 
the  average  results  shown  do  not  showcase  the  maximum  potential  of  MFC  in  the 
face  of  truly  extreme  wind  cases.  However,  because  the  average  bending  moment 
percent  reduction  in  all  cases  with  current  wind  information  is  well  over  50%,  these 
percent  reductions  can  only  increase  as  weather  conditions  become  less  stable.  This 
would  act  to  increase  the  launch  envelope  of  any  given  vehicle,  addressing  several  of 
the  motivating  concerns  described  in  Chapter  1. 

8.2  Focused  Optimized  MFC /Baseline  Performance 
Comparison 

The  two  main  purposes  of  this  comparison  set  are  to  comprehensively  showcase  the 
ability  that  bias  correction  signal  smoothing  has  to  significantly  decrease  gimbal  activ¬ 
ity  while  only  slightly  increasing  bending  moments  and  to  highlight  the  comparisons 
between  this  system  and  the  best  current  baseline  controller  available,  namely  the 
baseline  controller  using  a  15th  order  polynomial  fit  of  delayed  wind  profiles. 

As  shown  previously,  the  baseline  controller  performs  best  when  it  is  provided 
with  a  smoothed  version  of  the  actual  wind  profile.  This  controller /wind  processing 
approach  is  similar  to  current  leading  industry  lavmch  vehicle  control  concepts.  The 
smoothed  wind  profile  allows  the  gimbal  activity  to  remain  low,  while  still  providing 
enough  wind  information  to  the  controller  to  reasonably  address  a  and  Qa  concerns. 
This  combination  of  low  gimbal  activity  and  relatively  good  a  minimization  makes 
this  the  best  performing  basehne  combination  that  was  tested,  in  both  the  delay  and 
no  delay  cases.  When  paired  with  a  balloon  based  wind  measurement  system,  it 


141 


can  be  considered  equivalent  to  the  best  launch  vehicle  control  systems  in  existence 
today.  When  paired  with  a  LIDAR  wind  measurement  system,  it  would  hold  a  distinct 
advantage  over  other  launch  vehicle  control  systems  currently  in  use. 

In  an  effort  to  reduce  gimbal  activity  as  much  as  possible,  while  paying  the  lowest 
cost  in  bending  moment  performance,  both  the  ASAS  and  NSAS  controllers  were 
flown  against  the  complete  ETR  and  WTR  wind  profile  sets  while  smoothing  the  bias 
correction  signal  using  a  ±  4  second  moving  average.  The  results  follow,  in  figure  and 
table  form.  In  each  figure,  the  terms  Base  15  D  and  Base  15  ND  refer  to  the  baseline 
controller  using  a  15th  order  poljaromial  fit  of  the  delayed  and  actual  wind  profiles, 
respectively.  MFC  Inf  and  MFC  Inf  BC  refer  to  the  MFC  controllers  using  the  actual 
wind  profile  without  and  with  bias  correction  signal  smoothing.  The  tables  contain 
%  reductions  from  various  baseline  configurations  to  various  MFC  configurations.  As 
in  the  previous  set  of  tables,  a  positive  number  indicates  a  percent  reduction  and  a 
negative  number  indicates  a  percent  increase. 


Figure  8-11;  Feak  a  Response  to  Varying  Controller /Bias  Signal  Frocessing  Combi¬ 
nations 
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Figure  8-12:  Peak  Qa  Response  to  Varying  Controller/Bias  Signal  Processing  Com¬ 
binations 


Figure  8-13:  Peak  Bending  Moment  Response  to  Varying  Controller/ Bias  Signal  Pro¬ 
cessing  Combinations 
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Figure  8-14:  Peak  Gimbal  Angle  Deflection  Response  to  Varying  Controller/Bias 
Signal  Processing  Combinations 


Figure  8-15:  Total  Gimbal  Angle  Travel  Response  to  Varying  Controller/Bias  Signal 
Processing  Combinations 


a 

Qa 

B.M. 

G.A.D. 

T.G.A.T. 

AS  AS  ASAS  BP:  ETR 

-26.8706 

-26.8331 

-16.4674 

26.991 

14.3278 

NSAS  NSAS  BP:  ETR 

-15.3887 

-16.0626 

-5.1003 

32.2569 

27.859 

ASAS  ASAS  BP;  WTR 

-23.2472 

-23.6964 

-15.6449 

28.3301 

14.7906 

NSAS  NSAS  BP:  WTR 

-12.0968 

-13.2787 

-5.7497 

32.8116 

28.8208 

Table  8.2;  MFC  Metric  %  Reductions  Due  to  Bias  Correction  Signal  Smoothing 


These  results  show  that  bias  correction  signal  smoothing,  or  bias  processing  (BP) 
significantly  decreases  overall  gimbal  activity  in  both  MFC  controllers.  The  MFC 
controllers  using  smoothing  are  denoted  with  a  BP.  The  smoothing  produces  an  ap¬ 
proximate  30%  reduction  in  gimbal  angle  defiection  for  both  controllers.  Total  gimbal 
angle  travel  is  reduced  by  approximately  14%  in  the  AS  AS  controller  and  approxi¬ 
mately  28%  in  the  NS  AS  controller.  These  reductions  do  cause  a  shght  increase  in 
bending  moments,  with  larger  increases  in  both  a  and  Qa  metrics  as  well.  The  NS  AS 
controller  showed  the  smallest  increase  in  bending  moments,  at  approximately  5%, 
while  the  AS  AS  showed  an  approximate  increase  of  16%.  The  actual  results  of  the 
MFC  simulations  with  bias  signal  smoothing  is  included  in  Appendix  C.  As  stated 
previously,  steps  such  as  bias  signal  or  wind  processing,  may  or  may  not  be  necessary, 
depending  on  the  individual  vehicle  limitations.  While  gimbal  activity  is  not  a  factor 
for  the  K-1,  it  is  useful  to  show  that  the  tradeoff  can  be  beneficial  in  certain  cases. 
The  following  set  of  tables  show  that  this  smoothing  can  be  implemented  while  still 
producing  well  over  50%  reductions  in  a,  Qa  and  bending  moment  metrics. 

Tables  8.3  and  8.4  show  the  percent  reductions  which  occur  between  the  baseline, 
with  and  without  delays,  and  both  MFC  controllers,  with  and  without  bias  signal 
processing,  in  the  face  of  winds  from  both  test  ranges.  The  delay  and  no  delay  results 
will  be  dealt  with  separately. 
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a 

Qa 

B.M. 

G.A.D. 

T.G.A.T. 

Base  ND  ^  AS  AS:  ETR 

66.07 

65.13 

59.95 

-234.55 

-520.22 

Base  ND  NS  AS:  ETR 

-169.08 

-278.49 

Base  ND  AS  AS:  WTR 

-270.14 

-567.85 

Base  ND  ^  NS  AS:  WTR 

-195.50 

-304.82 

Base  D  ^  AS  AS:  ETR 

74.33 

-124.25 

-503.26 

Base  D  NSAS:  ETR 

74.02 

74.25 

77.71 

-80.36 

-268.14 

Base  D  ^  AS  AS:  WTR 

71.05 

71.10 

. - . - . - . 

66.96 

-175.77 

-561.82 

Base  D  ^  NSAS:  WTR 

70.82 

71.24 

75.44 

-120.15 

-301.16 

Table  8.3:  Metric  %  Reductions  Due  to  MFC  Controllers  without  Bias  Correction 
Signal  Processing 


a 

Qoc 

B.M. 

G.A.D. 

T.G.A.T. 

Base  ND  ^  AS  AS  BP:  ETR 

56.96 

55.77 

53.35 

-144.25 

-431.36 

Base  ND  ^  NSAS  BP:  ETR 

60.37 

59.32 

68.63 

-82.28 

-173.05 

Base  ND  ^  AS  AS  BP:  WTR 

54.24 

52.75 

50.24 

-165.28 

-469.07 

Base  ND  NSAS  BP:  WTR 

58.04 

56.93 

66.21 

-98.54 

-188.14 

Base  D  AS  AS  BP:  ETR 

67.44 

67.51 

65.17 

-63.72 

-416.82 

Base  D  ^  NSAS  BP:  ETR 

70.02 

70.12 

76.58 

-22.18 

-165.58 

Base  D  ^  ASAS  BP:  WTR 

64.32 

64.26 

61.74 

-97.64 

-463.93 

Base  D  NSAS  BP:  WTR 

67.28 

67.42 

74.02 

-47.92 

-185.54 

Table  8.4:  Metric  %  Reductions  Due  to  MFC  Conti'ollers  with  Bias  Correction  Signal 
Processing 

As  previously  stated,  the  baseline  controller  using  delayed  winds  represents  the 
best  controller/ wind  measurement  process  combination  currently  available  to  the 
space  launch  industry.  The  implementation  of  the  bias  signal  smoothing  malces  a 
very  small  impact  on  the  bending  moment  %  reduction  between  both  MFC  controllers 
and  the  baseline.  The  NS  AS  change  in  %  reductions  is  approximately  1%  while  the 
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AS  AS  change  is  approximately  5%.  With  the  bias  smoothing  in  place,  the  AS  AS  and 
NS  AS  controllers  still  produce  approximate  bending  moment  %  reductions  of  63% 
and  75%,  respectively,  over  the  basehne  controller  using  delayed  wind  data.  However, 
these  small  changes  in  the  bending  moment  metrics  produce  significant  changes  in 
the  %  that  the  gimbal  activity  metrics  change  between  the  basehne  and  the  MFC 
controllers.  The  implementation  of  the  bias  smoothing  causes  the  NS  AS  total  gim¬ 
bal  angle  travel  %  increase  over  the  baseline  to  change  from  approximately  285%  to 
175%,  a  greater  than  100  percentage  point  decrease.  Similarly,  an  approximate  65 
percentage  point  decrease  occurs  in  the  NS  AS  gimbal  angle  deflection  metric.  Re¬ 
ductions  such  as  this  exist  for  the  AS  AS  controller  as  well,  although  they  are  less 
prominent.  This  is  due  to  the  removal  of  direct  control  over  the  commanded  gimbal 
angle  from  the  MFC  controller  in  the  AS  AS  architecture. 

When  all  simulations  are  provided  with  current  wind  information,  the  results 
which  compare  the  baseline  and  MFC  controllers  highlight  the  advantages  that  MFC 
alone  has  to  offer  over  the  standard  reactive  control  system  approach.  When  the  bias 
signal  processing  is  implemented,  the  bending  moment  metric  %  reduction  experiences 
an  approximate  2  percentage  point  swing  in  the  NS  AS  controller  and  an  approximate 
6  percentage  point  swing  in  the  AS  AS  controller.  When  the  bias  smoothing  is  imple¬ 
mented,  the  average  bending  moment  reduction  in  the  AS  AS  and  NSAS  controller 
are  approximately  52%  and  67%,  respectively.  These  small  increases  in  bending  mo¬ 
ments  yield  significant  decreases  in  gimbal  activity.  As  a  result  of  implementing  bias 
smoothing,  the  NSAS  total  gimbal  angle  travel  %  increase  changes  from  approxi¬ 
mately  292%  to  approximately  180%,  a  112  percentage  point  decrease.  Similarly,  the 
NSAS  gimbal  angle  deflection  metric  %  increase  changes  from  approximately  182% 
to  approximately  91%,  a  91  percentage  point  decrease.  These  results  are  mirrored  in 
the  AS  AS  controller,  although  they  are  not  as  significant. 

This  research  shows  that  MFC  alone  can  reliably  yield  bending  moment  reductions 
of  approximately  58%  to  69%,  or  52%  to  67%  if  bias  signal  smoothing  is  implemented, 
depending  on  the  controller  layout.  Gimbal  activity  is  significant,  but  can  be  largely 
reduced  if  bias  signal  smoothing  is  implemented.  These  reductions  were  even  larger 
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when  the  MPC  controller’s  performance  was  compared  to  a  realistic  current  launch 
vehicle  control  approach.  Under  these  conditions,  the  MPC/LIDAR  combination 
could  reliably  yield  bending  moment  reductions  of  approximately  69%  to  76%,  or 
63%  to  75%  if  bias  signal  smoothing  is  implemented,  depending  on  the  controller 
layout. 


8.3  Wind  Measurement  Noise  Comparisons 

A  set  of  trials  was  completed  to  investigate  the  behavior  of  each  control  architecture  in 
the  face  of  increasing  noise  within  the  wind  measurement  process.  100  wind  profiles 
from  the  ETR  and  WTR  were  used  to  create  trajectorj^  sets  designed  for  winds 
containing  additional  wind  measurement  noise  signals  with  standard  deviations  of  0, 
2,  4,  8,  12,  16  and  20  ft/sec,  respective!}'.  Each  controller  was  then  flown  following 
the  trajectories  resulting  from  the  noisy  wind  profiles  while  only  being  affected  by  the 
actual  wind  profiles.  Current  wind  information  was  used  for  all  trials.  The  resulting 
performance  metrics  for  each  controller  are  shown  in  Figures  8-16  to  8-20. 


Wind  Uncertainty  Standard  Deviation  (ft/sec) 


Figure  8-16:  Controller  Peak  a  Response  to  Wind  Measurement  Noise 
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Figure  8-17:  Controller  Peak  Qa  Response  to  Wind  Measurement  Noise 


Figure  8-18:  Controller  Peak  Bending  Moment  Response  to  Wind  Measurement  Noise 
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Figure  8-19:  Controller  Peak  Gimbal  Angle  Response  to  Wind  Measurement  Noise 


Wind  Uncertainty  Standard  Deviation  (ft/sec) 


Figure  8-20:  Controller  Total  Gimbal  Angle  Travel  Response  to  Wind  Measurement 
Noise 


These  figures  show  the  general  wind  measurement  noise  sensitivity  of  each  con- 
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troller  type.  Because  the  baseline  controller  exhibits  such  poor  performance  in  no¬ 
noise  conditions,  it  performs  worse  for  all  levels  of  wind  measurement  noise.  However, 
relative  to  its  no-noise  performance,  the  baseline  controller’s  metrics  show  less  sen¬ 
sitivity  to  wind  noise  than  either  of  the  MFC  controllers.  Both  MFC  controllers 
exhibit  a,  Qa  and  bending  moment  metrics  which  are  much  smaller  than  the  equiv¬ 
alent  baseline  metrics  at  the  zero-noise  point.  But  as  noise  increases,  these  metrics 
increase  at  a  rate  that  is  higher  than  the  baseline  rate  of  increase  at  the  same  points 
along  the  wind  noise  curve.  However,  both  MFC  controllers  still  show  better  a,  Qa 
and  bending  moment  performance  than  the  baseline  controller,  even  at  the  highest 
wind  noise  conditions  tested.  In  fact,  the  MFC  controllers  show  nearly  the  same  or 
lower  a,  Qa  and  Bending  Moment  metric  values  at  wind  noise  standard  deviations 
of  20  ft/sec  as  the  baseline  controller  shows  at  wind  noise  standard  deviations  of  0 
ft /sec.  These  maximum  wind  noise  levels  far  exceed  any  projected  system  noise  con¬ 
ditions  which  might  exist  in  a  contemporary  wind  measurement  system.  Thus,  the 
MFC  controllers  are  shown  to  be  superior  to  the  baseline  controller  in  aU  conceivable 
wind  noise  situations. 
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Chapter  9 


Conclusions  and  Recommendations 

9.1  Conclusions 

As  described  in  Chapter  1,  this  thesis  proposes  a  two-pronged  approach  to  mini¬ 
mizing  uncertainty  and  lowering  launch  vehicle  metrics,  such  as  a,  Qa  and  bending 
moments.  The  first  prong  of  the  approach  is  Model  Predictive  Control,  which  lever¬ 
ages  the  accurate  and  timely  wind  disturbance  data  produced  by  the  second  prong 
of  the  approach,  a  LIDAR  wind  sensor.  As  shown  in  the  previous  chapter,  while 
each  of  these  can  provide  certain  benefits  by  themselves,  the  combination  of  the  two 
yields  benefits  far  superior  to  the  sum  of  their  individual  contributions.  Without 
the  accurate  and  timely  wind  disturbance  data  produced  by  the  LIDAR  wind  sensor, 
the  MPC  controller  cannot  accurately  predict  the  vehicle’s  performance,  and  loses 
much,  if  not  all,  of  its  advantage  over  the  baseline  controller.  Therefore,  these  two 
technologies  must  be  viewed  and  used  as  a  single  integrated  system,  replacing  both 
the  current  balloon-based  wind  measurement  system  and  the  current  launch  vehicle 
control  approach. 

With  the  recent  dramatic  advances  in  available  computational  power,  the  techni¬ 
cal  hurdles  impeding  the  realization  of  aerospace  apphcations  of  MPC  have  become 
less  and  less  significant.  Current  launch  vehicle  control  systems  possess  more  than 
enough  computational  power  to  implement  the  MPC  control  systems  described  in  this 
thesis,  assuming  the  matrix  computations  are  completed  before  launch.  The  simple 
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matrix  multiplications  and  additions  required  by  this  method  can  easily  be  completed 
during  flight  using  today’s  on-board  processors.  While  the  increase  in  computational 
requirements  is  significant  for  real  time  calculation  of  the  MPC  matrices,  near-term 
processors  exist  which  will  meet  and  surpass  these  needs  within  a  matter  of  years. 

By  describing  the  development  of  two  separate  and  unique  approaches  to  the 
MPC  control  problem,  represented  by  the  AS  AS  and  NS  AS  controllers,  this  thesis 
has  shown  that  MPC  is  a  flexible  approach  to  launch  vehicle  control.  The  AS  AS 
controller  has  demonstrated  that  MPC  can  be  adapted  as  an  outer  control  loop  around 
an  existing  stability  augmentation  system,  while  the  NS  AS  controller  has  showcased 
MPC’s  ability  to  handle  both  the  stability  augmentation  and  higher  level  guidance 
tasks,  simultaneously. 

The  AS  AS  controller  shows  great  promise  as  a  template  for  future  integrations  of 
MPC  control  into  current  launch  vehicle  control  systems.  As  shown  in  the  previous 
chapter,  the  performance  improvement  possibilities  with  this  t3q)e  of  control  system 
are  significant.  While  this  controller’s  high  level  of  gimbal  activity  may  be  an  issue 
for  some  launch  vehicles,  several  solutions  do  exist  for  lowering  this  activity,  such  as 
modifications  to  the  inner  SAS  or  some  form  of  a  feedback  signal  processing. 

The  NS  AS  controller  reduces  the  MPC  approach  to  its  most  simple,  direct  form. 
By  avoiding  the  “crutch”  of  an  inner  loop  SAS,  MPC  can  be  implemented  in  the 
most  efficient  manner  possible.  In  cases  where  the  MPC  matrix  calculations  are 
completed  before  launch,  this  system  is  feasible  using  today’s  technology.  Controllers 
which  feature  complete  real-time  MPC  computation  are  possible  using  near-term 
computational  resources.  Because  of  its  abilitj^  to  control  gimbal  activity,  maximum 
reductions  in  bending  moments  can  be  achieved  using  this  method. 

While  the  AS  AS  approach  appears  to  be  the  most  suitable  to  near  term  appli¬ 
cations  of  MPC,  the  NS  AS  controller  appears  to  hold  the  most  promise  in  the  long 
term,  in  the  area  of  vehicle  performance.  This  is  not  surprising,  as  this  configuration 
places  MPC  in  direct  control  of  the  input  command  to  the  plant.  In  future  launch 
vehicles,  with  undetermined  control  systems  and  sufficient  computational  power,  the 
direct  application  of  MPC,  in  the  form  of  the  NS  AS  controller,  seems  to  be  the  more 
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powerful,  efficient  and  productive  implementation  approach. 

Implementing  an  MPC/LIDAR  control  system  will  certainly  require  additional 
initial  investment  dollars  over  traditional  load  relief  and  wind  measurement  systems. 
However,  when  compared  to  the  benefits  this  system  can  deliver  to  the  typical  launch 
vehicle  operator  and  user,  these  costs  are  highly  justified  and  will  be  paid  back  to 
the  operator  in  a  relatively  short  time  span.  The  increased  launch  capability  and  risk 
reduction  inherent  to  the  MPC  approach  will  continue  to  yield  benefits  long  after 
the  additional  investment  has  been  recouped.  Significant  savings  will  also  be  realized 
by  replacing  the  current  manpower  intensive,  balloon-based  process  of  wind  measure¬ 
ment  with  the  more  accurate,  automated  and  timely  LIDAR-based  wind  measurement 
system. 

9.2  Recommendations  for  Future  Work 

As  explained  earlier,  both  MPC  controllers  exhibit  a  bias  during  flight  which  must  be 
corrected  in  order  to  achieve  the  maximum  load  relief  performance.  The  source  of  this 
error  is  uncertain  at  the  time  of  this  writing,  although  several  possible  explanations 
exist.  It  is  important  to  understand  that  the  error  may  merely  be  a  result  of  the 
unstable  plant  being  controlled  by  MPC.  This  error  does  not  exist  when  MPC  is  used 
to  control  a  stable  plant.  Ehminating  this  bias  would  be  useful,  since  it  would  allow 
the  simulation  to  be  flown  without  an  initiaUzation  cycle,  as  shown  in  this  thesis. 
The  initialization  cycle  is  an  acceptable  work-around  solution  to  this  drift  problem, 
but  certainly  not  the  most  elegant  or  efficient. 

Another  issue  which  should  be  investigated  is  a  more  realistic  development  of  the 
inner  SAS  in  the  ASAS  controller  configuration.  The  current  simulation  feeds  the 
actual  a  signal,  in  real  time,  to  the  a  SAS,  which  uses  it  to  calculate  the  next  gimbal 
angle  command  for  the  vehicle.  This  implicit  assmnption  of  accurate  a  knowledge 
may  not  be  completely  realistic.  Currently,  vehicle  a  knowledge  is  only  available  using 
estimation  techniques  of  uncertain  accuracy.  An  investigation  into  more  accurate 
modelling  of  these  methods  would  yield  important  performance  data  and  a  more 
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complete  understanding  of  the  performance  of  integrated  outer- loop  MPC  controllers 
with  current  inner-loop  stability  augmentations  systems. 

It  is  also  important  to  note  that  each  cm'rent  MPC  controller  parameter  was  se¬ 
lected  from  a  subset  of  the  complete  parameter  range  of  each  variable.  It  is  possible 
that  a  more  extensive  search  could  reveal  a  more  advantageous  parameter  set.  In  ad¬ 
dition,  each  MPC  controller  is  also  the  result  of  a  relatively  crude  selection  process, 
which  did  not  include  full  6DOF,  non-linear  models  or  full  ascent  wind  simulation. 
A  more  comprehensive  investigation  of  these  architectures  using  the  complete  simu¬ 
lations  could  yield  significantly  different  results  than  were  initially  recorded. 

The  controller  design  process  documented  in  this  thesis  operates  on  the  concept 
that  one  controller  parameter  set  can  be  found  which  can  be  used  for  virtually  any 
possible  wind  profile.  While  this  was  demonstrated  to  be  true  through  the  use  of  over 
3000  actual  wind  profiles  occurring  over  nearly  four  decades,  the  subject  of  designing 
the  controller  for  an  individual  wind  profile  was  not  investigated.  By  definition,  the 
final  controller  parameter  set  was  the  average  optimum  result  for  the  127  profile  design 
wind  set.  Even  larger  performance  gains  are  possible  if  this  modelling  was  done  on 
an  individual  basis.  This  process  could  conceivably  be  automated  and  condensed  into 
a  process  of  only  a  few  minutes  or  less,  enabling  maximum  controller  performance  in 
the  face  of  any  given  wind  profile. 


156 


Appendix  A 


Notation  Convention  and  Variable 
Decl£iration 

A.l  Notation 

Throughout  this  thesis,  very  specific  notations  conventions  are  utilized.  These  con¬ 
ventions  allow  the  precise  and  succinct  labelling  of  variables  and  constants. 

1.  Transformation  matrices  are  referred  to  using  the  following  convention: 


TgJ  where:  Sf  =  Source  frame 

Df  =  Destination  frame 

An  example  of  this  naming  convention  is  shown  below: 

Tg  =  The  transformation  matrix  from  the  earth  frame  to  the  inertial  frame. 

In  addition,  individual  positions  in  a  transformation  matrix  are  referred  to  using 
their  row  and  column  numbers,  enclosed  in  brackets,  T[row][column]. 
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An  example  of  this  naming  convention  is  shown  below: 

r[2][3]  =  The  term  located  in  the  second  row  and  in  the  third  column  of  the 
transformation  matrix,  T. 

2.  Vectors  are  referred  to  using  the  following  naming  convention: 


where:  V  — Vector  type 

Vs  =  Vector  source 
Frwt  =  Frame  with  respect  to 
Fex  =  Frame  expressed  in 

An  example  of  this  naming  convention  is  shown  below: 

V  =  The  acceleration  of  the  body,  with  respect  to  the  earth  frame,  expressed 
in  the  inertial  frame. 


A. 2  Variable  Declaration 

Throughout  this  thesis,  many  variable  designations  have  been  used.  The  description 
and  units  of  these  are  included  below.  Note  that  although  the  standard  angular 
units  in  the  simulations  are  radians,  degrees  were  used  when  analyzing  the  data  for 
simplicity. 

Variable  Designation  Variable  Description  Variable  Units 


Vector  Type 

Ij  Acceleration 

Td  Angular  Acceleration 

Zj  Angular  Velocity 


ft/sec^ 
rad/sec^ 
ft /sec 
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F 

Force 

Ibf 

M 

Moment 

ft-lb 

'P 

Position 

ft 

V,v 

Velocity 

ft /sec 

Vector  Source 

Aref 

(see  below) 

- 

a 

Air 

- 

b 

Body 

- 

e 

Earth  Frame 

- 

eng 

Engine 

g 

Gravity 

Ig 

Local  Geographic  Frame 

Frames 

Aref 

Body  Frame  centered  at  Aero  Reference  Point 

_ 

B,  b 

Body  Frame 

- 

I,  i,  ECI 

Inertial  Frame 

E,  e,  ECR 

Earth  Frame 

LG,  Ig,  NED 

Local  Geographic  Frame 

Aerodynamic  Terms 

a 

Angle  of  Attack 

rad 

Total  Angle  of  Attack 

rad 

P 

Sideslip  Angle 

rad 

Ca 

Axial  Force  Coefficient 

- 

Cy 

Lateral  Force  Coefficient 

- 
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Normal  Force  Coefficient 

- 

C'„^ 

Pitching  Moment  Coefficient 

- 

Cl 

Rolling  Moment  Coefficient 

- 

Cn 

Yawing  Moment  Coefficient 

- 

L 

Reference  Length 

ft 

P 

Roll  Rate 

rad/sec 

(t> 

Euler  Roll  Angle 

rad 

4>* 

Aerodynamic  Roll  Angle 

rad 

psi 

Euler  Yaw  Angle 

rad 

Q 

Pitch  Rate 

rad/sec 

Q,  Qbar 

Dynamic  Pressure 

Ibf/ft^  (psf) 

r 

Yaw  Rate 

rad/sec 

P 

Air  Density 

slug/ft^ 

S 

Reference  Area 

ft2 

SoS 

Speed  of  Sound 

ft  /  sec 

e 

Euler  Pitch  Angle 

rad 

Round  Earth  Terms 

GMST 

Greenwich  Mean  Solar  Time 

sec 

k 

Earth  Flattening  Constant 

- 

Alt 

Altitude 

ft 

Lat 

Latitude 

rad 

Lon 

Longitude 

rad 

re 

Local  Earth  Radius 

ft 

reO 

Equatorial  Earth  Radius 

ft 

Engine  Terms 

Ae  Nozzle  Exit  Area  ft^ 
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Commanded  Gimbal  Angle 

rad 

Sp 

Pitch  Gimbal  Angle 

rad 

5y 

Yaw  Gimbal  Angle 

rad 

Prp 

Total  Engine  Thrust 

Ibf 

Ispvac 

Vacuum  Specific  Impulse 

sec 

m 

Mass  Flow  Rate 

slugs/sec 

TS 

Throttle  Setting 

- 

Tvac 

Vacuum  Thrust 

Ibf 

Mass  Property  Terms 

I XXI  lyyi  ^zz 

Body  Moments  of  Inertia 

slug-ft^ 

^xyt  ^yzi  Ixz 

Body  Cross-Products  of  Inertia 

Slug-ft2 

CG 

Center  of  Gravity 

- 

CGx 

CG  Location  in  the  Body  X  Axis  Direction 

ft 

General 

ASAS 

MPC  Controller  w/  an  o  SAS 

- 

Base 

Baseline  PID  Controller 

- 

ETR 

Eastern  Test  Range,  Kennedy  Space  Center,  FL 

- 

KAC 

Kistler  Aerospace  Corporation 

- 

LIDAR 

Light  Detection  and  Ranging 

- 

LQR 

Linear  Quadratic  Regulator 

- 

LTI 

Linear  Time  Invariant 

- 

MPC 

Model  Predictive  Control 

- 

NSAS 

MPC  Controller  w/o  a  SAS 

- 

PID 

Proportional-plus-Integrator-plus-Derivative 

- 

SAS 

Stabihty  Augmentation  System 

- 

TWD 

Tail-Wags-Dog 

- 
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WTR 

6D0F 


Western  Test  Range,  Vandenburg  AFB,  CA 
6  Degree  of  Freedom 


Appendix  B 


Simulation  S-Functions 


This  appendix  contains  a  short  summary  of  the  equations  executed  by  each  S- Function. 
All  necessary  matrix  transpose  operations  can  be  assumed  and  will  not  be  specifically 
mentioned.  All  matrix  equations  will  use  the  T[row]  [column]  convention  illustrated 
in  Appendix  A.  All  symbols  and  vector  conventions  are  also  explained  in  Appendix 
A.  All  angles  are  modified  to  fall  below  27r. 

accel 

m 

^(-^a+englB  “  ^  (B-2) 

eulerRates 


Ip  =  [§sin(0)  +  rcos(^)]/cos(^)  (B.3) 

0  =  g' cos(0)  —  r  sin(0)  (B.4) 

0  =  p  +  gtan(^)  sin(0)  +  rtan(0)  cos((/))  (B.5) 
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envGMST 


GMSTZ  =  GMSTZ{t  =  0)  +  a;e  *  Atim.e  (B.6) 

This  angle  is  used  to  perform  a  simple  rotation  about  the  Z  axis  to  produce  the 
current  T^. 

xformEuler 


2^^[l][l]  =  cos{6)  cos{ip) 

Tf[l][2]  =  cos(0)  sin('0) 

I-'llIP)  =  -sin(#) 

Tf  [2]  [1]  =  sin(<^)  sin(0)  cos(^/')  —  cos((j!>)  sin('0) 

jf[2][2]  =  sin(0)  sin(0)  sin(^>)  +  cos((5f))  cos('?/-’)  (B.7) 

2]^[2][3]  =  sin(0)cos(^) 

Tf  [3]  [1]  =  cos(<^)  sin(0)  cos(V-’)  +  sm{4>)  sin(V’) 

Tf[3][2]  =  cos(<3!))  sin(0)  sin(V’)  —  sin(^)  cos(V.!) 

Tf[3][3]  =  cos(<?>)  cos(^) 


t:  =  TtTl 


(B.8) 
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velocityOmega 


i?p[iiiii 

rjlilH 

rjlilH 

r,'[2][i] 

Ttg[m 

rjisiii] 

Tt,im 

I’lpPliai 


—  cos{Lon)  sm{Lat) 

—  sin(Lon) 

—  cos  (Lon)  cos  (Lot) 

—  sin(Lon)  sin(Lai) 
cos(Lon) 

—  sin(Lon)  cos(Lat) 
cos(Lat) 

0 

—  sin(Lat) 


(B.9) 


>7~i6  _  rpbrpe 

■^Ig  - 


(B.IO) 


elB 


■TI-^Wb 

(B.ll) 

[0  0  tOe] 

(B.12) 

=  Tf’uj^  1 

-^i  e  1/ 

(B.13) 

1  1 

(B.14) 

+  ^e|/  'Pb\i 

(B.15) 
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Lat  = 


h>*  ^6U(-) - ^ - ;;; - ^ 

kh'l  +  'pi\j{zY 


Lon  = 


'pI\,(^)  *  ^lUfa)  ~  ♦  ^UeW 


(B.16) 

(B.17) 

(B.18) 


=  Lai  *  sin(Lon) 

tafg|^[2]  =  —Lat  *  cos{Lon.)  (B.19) 

=  Lon 


— >/(7  — 

^  b  ^9  B 


Ve  j  =  iPb 


If*  = 

^  e  B  ^  e  j 


Vl  V  \  B~  '^ej 


(B.20) 

(B.21) 

(B.22) 

(B.23) 

(B.24) 

(B.25) 
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gravityShepperd 
This  model  produces  ~v 


aeroPeirameters 


V 


-^Ig^  a  lG 


B 


(B.26) 

(B.27) 

(B.28) 

(B.29) 


^  following  the  procedures  outlined  in  [8]  and  [9]. 

^  sIb  ~  oIe 

weight  =  m\v  g\ 


(B.30) 


(B.31) 


(B.32) 


V  =  \irtU  (B-33) 

l3=^sm-\irt\s{2)/V)  (B.34) 

a^  =  cos-\irt\^{l)/V)  (B.35) 
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X 


Figure  B-1:  Aerodynamic  Angles 

Q  =  tan“^Cu^|g(3)/‘u^|^(l)) 

(B.36) 

(j)*  -  tan-i(V^|^(2)/V^|^(3)) 

(B.37) 

mach  —  V/SoS 

(B.38) 

Q  =  l/2p  1/2 

(B.39) 

atmosphereU  S  76 

This  is  a  1976  U.S.  standard  atmospheric  model,  as  detailed  in  [6]. 
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xyz21at_lon_alt 


re 


1 


/ 


\ 


reO 


1  + 


pl\r*Pl\ 


k-1 


(B.40) 


engine 


M=  I'pJ  „|  -  re 


(B.41) 


Lon  —  arctan 


■p  6  e(^) 


(B.42) 


/ 


Lot  =  arctan 


^6  b(^) 


\ 


+  "p  6|£;(p)^)  kcos{Lat) 


Altil  -  k)  (B.43) 


m  = 


TS  *  Tvac 

Ispvac*  [Vc 


(B.44) 


Ft  =  TS  *  Tvac  —  Air  Pressure  *  Ae 


(B.45) 


^  eng\^[^] 

=  cos(5p)  cos{5y)FT 

F  eng  \  q  [v] 

=  —  cos{5p)  sin{5y)FT 

(B.46) 

=  sm((5p)  cos{5y)FT 

Aleng  ^  —  Peng  b  ^  ^ eng  b  (B.47) 
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aeroKlBoost 


This  function  uses  linear  interpolation  to  select  normal,  axial  and  lateral  force  co¬ 
efficients,  and  roll,  pitch  and  yaw  moment  coefficients  using  mach  number  and  q^. 
These  coefficients  are  denoted  as  Cn,  Ca,  Cy,  Cl,  Cm  and  Cn  respectively. 


Fn 

=  QSCn 

(B.48) 

Fa 

=  QSCa 

(B.49) 

Fl 

=  QSCy 

(B.50) 

Mr 

=  QSLCl 

(B.51) 

Mp 

=  QSLCm 

(B.52) 

My 

=  QSLCn 

(B.53) 

'Fa\ 

=  Fa 

Fa\ 

Ib[^] 

=  -Fl  cos{(l)*)  +  Fn  sm((f)*) 

(B.54) 

'Fa\ 

IbN 

=  Fl  sin{(i)*)  +  Fn  cos{(p*) 

=  Mr 

M  1  g  [p]  =  Mp  cos{(j)* )  + My  sin(0* ) 

(B.55) 

Ma^'^^Jy]  =  -  Mp  +  My  cos{(f)*) 

+  'PAre,\B  ^  (B-5S) 
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Appendix  C 


Simulation  Results 


The  results  shown  in  this  appendix  are  from  the  complete  wind  profile  set,  1,141  wind 
profiles  from  the  WTR  and  1,927  wind  profiles  from  the  ETR. 


C.l  Full  Wind  Set  Comparison  Results 

Results  from  ETR  Winds  without  Delay 


Base  6 

Base  15 

Base  Inf 

AS  AS  Inf 

NS  AS  Inf 

a  Avg. 

1.268 

1.0109 

1.0657 

0.34298 

0.34721 

a  St. Dev. 

0.41653 

0.28236 

0.28821 

0.083768 

0.085492 

Qoi  Avg. 

544.0496 

430.2306 

460.4586 

150.0318 

150.7934 

Qa  St. Dev. 

182.8657 

124.8899 

132.2233 

37.7273 

38.0168 

B.M.  Avg. 

581133.8751 

464495.7011 

568080.7079 

186039.58 

138646.6724 

B.M.  St.Dev. 

191880.8066 

131364.4525 

152070.0949 

45653.7213 

33862.486 

G.A.D.  Avg. 

0.16908 

0.14927 

0.5476 

0.49939 

0.40166 

G.A.D.  St.Dev. 

0.066439 

0.068349 

0.13759 

0.12445 

0.10552 

T.G.A.T.  Avg. 

4.6018 

4.6044 

21.2122 

28.5573 

17.4272 

T.G.A.T.  St.Dev. 

0.75429 

0.79842 

3.3997 

5.9521 

3.1971 
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Results  from  WTR  Winds  without  Delay 


Base  6 

Base  15 

Base  Inf 

AS  AS  Inf 

NS  AS  Inf 

Q  Avg. 

1.3115 

1.0503 

1.1014 

0.38996 

0.39315 

Q  St. Dev. 

0.51623 

0.35849 

0.3748 

0.17667 

0.17374 

Qa  Avg. 

570.4877 

460.4039 

495.2585 

175.8712 

175.0687 

Qa  St. Dev. 

241.2194 

173.9208 

185.9545 

74.6883 

73.0384 

B.M.  Avg. 

601636.9623 

484189.9257 

592250.4283 

208352.6563 

154689.7569 

B.M.  St.Dev. 

240990.6595 

171837.667 

199489.1667 

95371.7645 

73042.8821 

G.A.D.  Avg. 

0.1767 

0.14942 

0.58775 

0.55308 

0.44154 

G.A.D.  St.Dev. 

0.088578 

0.062783 

0.17337 

0.18183 

0.14927 

T.G.A.T.  Avg. 

4.4559 

4.4122 

22.0244 

29.4669 

17.8614 

T.G.A.T.  St.Dev. 

0.78769 

0.80865 

4.1171 

7.25 

3.8455 

Results  from  ETR  Winds  with  Delay 


Base  6 

Base  15 

Base  Inf 

AS  AS  Inf 

NS  AS  Inf 

a.  Avg. 

1.4168 

1.3363 

1.4656 

0.85738 

1.3106 

a  St.Dev. 

0.59751 

0.63388 

0.71057 

0.43991 

0.86828 

Qa  Avg. 

623.5383 

585.7086 

647.8754 

377.8147 

586.8081 

Qa  St.Dev. 

279.5638 

282.389 

317.9526 

186.7863 

383.8624 

B.M.  Avg. 

651627.9908 

622113.4911 

732958.6135 

468970.4965 

640633.4827 

B.M.  St.Dev. 

280245.4606 

291825.462 

333230.6556 

209664.483 

410049.2085 

G.A.D.  Avg. 

0.19463 

0.2227 

0.5944 

0.66837 

0.49532 

G.A.D.  St.Dev. 

0.10313 

0.14074 

0.19454 

0.24516 

0.20646 

T.G.A.T.  Avg. 

4.5435 

4.7.339 

22.7386 

37.1063 

19.289 

T.G.A.T.  St.Dev. 

0.8316 

0.97989 

4.2621 

8.0159 

3.9941 
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Results  from  WTR  Winds  with  Delay 


Base  6 

Base  15 

Base  Inf 

AS  AS  Inf 

NS  AS  Inf 

a  Avg. 

1.4447 

1.3471 

1.5102 

0.88017 

1.3462 

a  St. Dev. 

0.68639 

0.67906 

0.74221 

0.46195 

0.91945 

Qa  Avg. 

645.9825 

608.6365 

680.7474 

394.4207 

612.053 

Qa  St. Dev. 

318.9512 

309.9272 

347.4357 

201.6737 

412.9823 

B.M.  Avg. 

667723.2229 

629766.3976 

758800.298 

485601.4973 

661437.2946 

B.M.  St. Dev. 

319384.6031 

318289.0028 

356519.2014 

221628.7625 

438570.0142 

G.A.D.  Avg. 

0.20447 

0.20056 

0.62002 

0.6856 

0.51599 

G.A.D.  St.Dev. 

0.11554 

0.12238 

0.21711 

0.24581 

0.21478 

T.G.A.T.  Avg. 

4.4089 

4.4524 

22.1463 

35.1124 

18.6623 

T.G.A.T.  St.Dev. 

0.84254 

0.89813 

4.4773 

7.9057 

4.1409 
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C.2  Focused  Optimized  MFC /Baseline  Performance 
Comparison  Results 

Results  of  MFC  Controllers  with  Bias  Signal  Processing 


ASAS-.ETR 

NS  AS  :ETR 

ASAS-.WTR 

NSAS-.WTR 

a  Avg. 

0.43514 

0.40064 

0.48062 

0.44071 

a  St. Dev. 

0.10277 

0.096555 

0.18169 

0.17674 

Qa  Avg. 

190.2901 

175.0148 

217.5464 

198.3155 

Qa  St. Dev. 

47.0706 

43.8483 

80.0166 

76.3696 

B.M.  Avg. 

216675.4129 

145718.1052 

240949.144 

163583.9486 

B.M.  St.Dev. 

53240.2255 

36175.5857 

98192.2195 

74291.8856 

G.A.D.  Avg. 

0.3646 

0.2721 

0.39639 

0.29667 

G.A.D.  St.Dev. 

0.088407 

0.07436 

0.14833 

0.10215 

T.G.A.T.  Avg. 

24.4657 

12.5722 

25.1086 

12.7136 

T.G.A.T.  St.Dev. 

5.6506 

2.3763 

6.7018 

2.799 

Appendix  D 


Simulation  Parameters 


The  variables  contained  in  the  simulation  initialization  file  are  shown  in  this  appendix. 
The  variables  used  to  initialize  the  mass  properties  block  are  listed  in  Chapter  2. 

a  SAS  Controller  Gains 


Variable  Name  Variable  Value  Variable  Units 


Ka 

1.4 

- 

Kd 

1.2 

sec 

Kp 

0.7 

- 

Ks 

0.1 

_ 

Initial  State  Values 


Variable  Name 

Variable  Value 

Variable  Units 

Initial  Position 

[0  0  0] 

ft 

Initial  Velocity 

[0  0  0] 

ft/sec 

Initial  Angular  Position 

[0  0  0] 

rad 

Initial  Angular  Velocity 

[0  0  0] 

rad/sec 

Initial  Mass 

26243 

slugs 
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Engine  Parameters 


Variable  Name 

Variable  Value 

Variable  Units 

Engine  Position 

[34.2667  0  0] 

ft 

Vacuum  Thrust 

378000 

Ibf 

Vacuum  Isp 

331.3 

sec 

Exit  Area 

18.479167 

ft2 

Aerodynamic  Parameters 

Variable  Name 

Variable  Value 

Variable  Units 

Aero  Reference  Point 

[83.333  0  0] 

ft 

Reference  Area 

201.0 

ft^ 

Reference  Length 

16.0 

ft 

Actuator  Dynamics  and  Atmospheric  Variables 

Variable  Name 

Variable  Value 

Variable  Units 

25.132 

rad/sec 

C 

0.707 

- 

Earth  Radius 

20925741 

ft 

Gravity  Inverse 

0.03108 

sec^/ft 
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Appendix  E 


SIMULINK  Non-Linear  6DOF 
Implementation 


E.l  SIMULINK  Model  Structure 

The  SIMULINK  operating  environment  allows  the  user  to  create  complex  simula¬ 
tions  with  any  number  of  subsystems  nested  in  multiple  levels.  This  enables  the 
implementation  of  highly  complicated  system  simulations  while  still  maintaining  a 
simple,  easy  to  understand  interface.  The  following  sections  contain  pictures  of  the 
significant  systems  and  subsystems  contained  within  the  NS  AS,  AS  AS  and  Baseline 
models.  Because  the  simulations  are  identical  except  for  the  pitch  controller  in  each 
model,  only  one  set  of  of  pictures  will  be  included  for  the  remainder  of  the  Rocket 
subsystems,  which  are  located  in  the  NS  AS  Model  section. 

The  Scopes  subsystem  exists  in  all  three  simulations,  and  contains  SIMULINK 
variable  plotting  blocks  which  can  be  used  to  observe  different  signals  during  or  after 
simulation.  The  three  outputs  from  the  Scopes  subsystem  are  the  bending  moments, 
calculated  at  stations  located  at  various  points  on  the  vehicle.  The  station  used  for 
all  bending  moment  readings  in  this  thesis  was  the  furthest  aft  station.  Station  1000, 
located  approximately  at  the  CG. 
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E.2  NS  AS  Model 


The  Wind  XYZ  subsystem  contains  the  actual  wind  disturbance  profile  which  affects 
the  vehicle  during  flight  and  is  common  to  all  three  simulations.  The  Actuator  Dy¬ 
namics  block  is  a  simple,  second  order  transfer  function  using  the  previously  deflned 
values  of  and  (^. 


Figure  E-2:  NS  AS  Model:  Rocket 


The  Rocket  subsystem  contains  two  subsystems,  Forces  and  Moments  and  Vehicle 
Motion,  as  well  as  a  flat-earth  version  of  the  S-function,  velocityOmega_SMLK.  bf- 
fCGPos  is  the  position  of  the  CG  in  the  body  frame.  altAhead  is  a  vector  of  projected 
vehicle  altitudes.  The  vector  contains  p  steps,  with  a  step  size  equal  to  the  prediction 
step  size.  Outputs  are  shown  as  ovals.  If  these  exist  in  the  top  level  of  a  model,  the 
signals  are  passed  out  to  the  Matlab  workspace.  Otherwise,  they  pass  signals  up  to 
the  next  subsystem  level. 
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Figure  E-3:  NSAS  Model:  Pitch  Controller 


The  NSAS  pitch  controller  contains  both  the  MPC  and  PID  controllers.  The 
switch  transfers  control  from  the  PID  controller  to  the  hIPC  controller  at  a  time 
corresponding  to  passage  of  the  appropriate  Q  level.  Because  the  research  is  only 
conducted  in  the  pitch  plane,  which  is  also  the  inertial  Z  plane,  the  wind  disturbance 
profiles  is  contained  in  a  look-up  table  titled  Z  wind  Dist.  This  table  produces  a 
vector  of  disturbances  from  each  altAhead  signal,  which  is  then  loaded  into  the  MPC 
controller. 
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Figure  E-4:  NS  AS  Model:  Forces  and  Moments 

The  Forces  and  Moments  subsystem  contains  all  the  S-functions  required  to  cal¬ 
culate  the  forces  and  moments  due  to  engine  forces,  aerodynamic  forces  and  gravity. 
These  S-functions  include  atmosphereUS76_SMLK,  engine_SMLK,  aeroKlBoost_SMLK, 
gravityShepperd-SMLK  and  AeroParameters_SMLK.  The  total  forces,  neglecting  the 
effects  of  gravity  are  also  foimd,  for  use  in  calculating  the  lateral  acceleration  of  the  ve¬ 
hicle.  Gains  are  modelled  as  triangles  and  are  used  most  frequently  in  this  subsystem 
to  change  units  from  radians  to  degrees,  using  the  variable  r2d. 
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Figure  E-5:  NS  AS  Model:  Vehicle  Motion 


The  Vehicle  Motion  subsystem  contains  the  S-functions  required  to  propagate  the 
angular  position  of  the  vehicle  through  time,  eulerRate_SMLK  and  a  flat  earth  version 
of  xformEuler_SMLK.  It  also  contains  the  S-function  which  produces  the  altAhead 
vector  of  projected  altitudes,  altAhead_SMLK. 
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Figure  E-6:  NSAS  Model:  Motion 


The  Motion  subsystem  contains  the  S-function,  acceLSMLK,  which  calculates 
hnear  and  rotational  accelerations.  The  integrator  blocks  used  to  calculate  the  hnear 
velocity  and  position  and  the  rotational  velocity  are  also  visible.  In  addition,  the 
lateral  acceleration  is  calculated  using  the  no-gravity  forces  calculated  in  the  Forces 
and  Moments  subsystem.  The  Mass  Properties  subsystem  calculates  the  inertia  tensor 
and  CG  position  as  they  change  over  time. 
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E.3  AS  AS  Model 


Figure  E-7;  ^5^5  Model;  Top  Level 


The  AS  AS  simulation  top  level  looks  verj^  similar  to  the  NS  AS  simulation  top  level 
with  the  exception  of  the  additional  integrator  state  which  is  passed  into  the  pitch 
controller. 
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Figure  E-8:  ASAS  Model:  Pitch  Controller 


The  ASAS  pitch  controller  subsystem  is  very  similar  to  the  NS  AS  pitch  con¬ 
troller,  with  the  addition  of  the  a  SAS  between  the  MPC  controller  and  the  output 
switch.  The  block  labelled  no  int  WU  is  used  to  eliminate  integrator  windup  during 
the  time  when  the  MPC  controller  is  not  operating. 
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E.4  Baseline  Model 


Figure  E-9:  Baseline  Model:  Top  Level 

The  Baseline  simulation  shows  the  common  subsj^stems  alreadj^  discussed  in  previous 
sections.  The  pitch  and  yaw  control  are  both  contained  in  the  subsystem  Classical 
G&C. 
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Cl>— ► 


Yaw  Euler  Angle 


Yaw  delta  Command 


Yaw  Euler  Angle 


Yaw  GImble  Command 


Yaw  Controller 


Figure  E-10:  Baseline  Model:  Classical  G&C 


The  Baseline  pitch  controller  is  shown  above.  The  profile  loaded  into  the  look-up 
table  Theta  Traj  is  the  pitch  trajectory  calculated  by  the  trajectory  designer.  The 
pitch  controller  follows  this  commanded  trajectory.  This  controller  implements  the 
gain  scheduling  described  in  chapter  6. 
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